Improved microscopic-macroscopic approach incorporating the 

effects of continuum states 



Naoki Tajima 
Department of Applied Physics, University of Fukui, 
3-9-1 Bunkyo, Fukui 910-8507, Japan 

Yoshifumi R. Shimizu 

Department of Physics, Graduate School of Science, 
Kyushu University, Fukuoka 812-8581, Japan 

Satoshi Takahara 

Kyorin University, School of Medicine, Mitaka, Tokyo 181-8611, Japan 

(Dated: June 14, 2010) 

Abstract 

The Woods-Saxon-Strutinsky method (the microscopic-macroscopic method) combined with 
Kruppa's prescription for positive energy levels, which is necessary to treat neutron rich nuclei, 
is studied to clarify the reason for its success and to propose improvements for its shortcomings. 
The reason why the plateau condition is met for the Nilsson model but not for the Woods- 
Saxon model is understood in a new interpretation of the Strutinsky smoothing procedure as a 
low-pass filter. Essential features of Kruppa's level density is extracted in terms of the Thomas- 
Fermi approximation modified to describe spectra obtained from diagonalization in truncated 
oscillator bases. A method is proposed which weakens the dependence on the smoothing width 
by applying the Strutinsky smoothing only to the deviations from a reference level density. The 
BCS equations are modified for the Kruppa's spectrum, which is necessary to treat the pairing 
correlation properly in the presence of continuum. The potential depth is adjusted for the 
consistency between the microscopic and macroscopic Fermi energies. It is shown, with these 
improvements, that the microscopic-macroscopic method is now capable to reliably calculate 
binding energies of nuclei far from stability. 
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I. INTRODUCTION 



Understanding the properties of unstable nuclei is one of the most interesting subjects 
of nuclear physics It is also important for astrophysics; for example, determination of 
the precise position of neutron drip line is crucial for the r-process nucleosynthesis A 
characteristic feature of unstable nuclei, among others, is the weak binding of nucleons, so 
that the proper treatment of continuum (scattering) states is very important for the two 
basic ingredients of the nuclear structure, the shell effect and the pairing correlation j^|. 
The most popular method of recent years to treat this problem is the selfconsistent mean 
field theory, especially the Hartree-Fock-Bogoliubov (HFB) theory fl, with suitably cho- 
sen (density dependent) zero- or finite-range effective interactions [5]. Such selfconsistent 
mean field models can reproduce the very basic quantities like the nuclear mass rather 



well jgl, and can be used to investigate the detailed deformation properties of nucleus. 
On the other hand, a non-selfconsistent semi-phenomenological method of the Strutinsky 
shell correction approach jz iO], or often called the microscopic-macroscopic method, has 
been used for more than forty years in order to calculate nuclear masses, deformations 
and fission paths. In such an approach, the part of binding energy smoothly varying as a 
function of nucleon (proton and neutron) number is represented by the liquid-drop or the 
droplet model with parameters adjusted to reproduce the experimental binding energy, on 
top of which is added the rapidly varying shell energy correction evaluated by assuming 
some non-selfconsistent single-particle potential. 

It is known that there is a close relationship between the two, the selfconsistent mean 



field and the shell correction approaches III, |12J], but the actual calculational procedures 



differ considerably and their own merits are quite different. The number of adjustable 
parameters is generally fewer and the range of applicability is believed to be wider in 
the selfconsistent mean field models, whereas the shell correction approach requires much 
less computational power. Thanks to recent progress of computer power, the root mean 
square deviation between the calculated and experimental masses, as an example, in some 



of the selfconsistent mean field models 
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14j is approaching the same level of accuracy 



as that in the state of the art model of the shell correction approach [15j| (or even better). 
However, its ease of computation and its flexibility of choosing the single-particle potential 
are still great merits of the shell correction approach. For example, the various effects of 



the single-particle orbits can be more directly studied in the shell correction approach. 
On the contrary, in selfconsistent mean field models, a clear-cut picture is sometimes lost 
due to the complicated selfconsistency between the nuclear mean field and the effective 
interaction. 

Although the qualities of the mass fit are similar in the two approaches in stable nuclei, 
they often give quite different predictions for very heavy nuclei and unstable nuclei near 
the neutron drip line, where no experimental date are available jf], [l^, [ijj]. It should be 
noticed that the shell correction approach has several difficulties for the calculation of 
unstable nuclei, which are mainly related to the problem of unbound (continuum) states 
characteristic to weakly bound systems. The difficulties were carefully examined one by 
one in Ref. 16]. 

The first and the most crucial difficulty is that the shell correction energy cannot 
be unambiguously determined for the single-particle potential with finite depth, which 
is indispensable for describing weakly bound states. The energy of shell correction is 
defined as the difference between the sum of single-particle energies up to the Fermi level 
and its smoothed part. The conventional way of the smoothing procedure utilizes the 
energy averaging of the single-particle level density over the interval 7 of the typical shell 
spacing, 7 « %uj = Al/A 1 ^ MeV, where A is the mass number. If the absolute energy 
of the Fermi level is smaller than the shell spacing, the averaging inevitably involves the 
unbound states. The continuum single-particle levels are usually discretized by using, e.g., 
the harmonic oscillator basis expansion, but blind inclusion of them leads to divergent 

nn 

results as the basis is enlarged even in stable nuclei [16|, 1171] ; this is simply because the 
level density of continuum states itself is a divergent quantity. It is proposed that the 
level density above the threshold should be replaced [18| with the so-called continuum 



level density 19|, |20(, and the resultant smoothed energy is shown to be convergent 21]. 



The evaluation of the continuum level density requires the energy derivative of the 
phase shift, or of the scattering matrix in general, so that the calculation is cumber- 
some for spherical nuclei and difficult for deformed nuclei. A breakthrough was given by 
Kruppa [22], who proposed a powerful practical prescription to calculate the continuum 
level density by using the fact that it is written as the difference between the level den- 



sities with and without the mean field potential 
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241 ] . However, the problem remains; 
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the so-called plateau condition [25], which guarantees that the shell correction energy is 



independent of the smoothing procedure, is not well satisfied generally [26|, |27| . We rein- 
vestigate the meaning of energy smoothing procedure and consider a remedy to recover 
the plateau condition as much as possible employing the Kruppa's prescription. 

It is worth mentioning that there are different methods to calculate the smoothed 



part. One is to make averaging with respect to the partic 



e number, not to the single- 



particle energy, only by employing the bound states [28|, |29|. However, the resultant 
smoothed energy depends sensitively on how to perform the averaging for nuclei near 
the drip line 3CH33| . where there is no unoccupied bound states and thus one has to 
tackle a difficult task to estimate an average value at a point (a particle number) using 
data points only on one side of that point (at smaller particle numbers). It is also a 
problem that the smoothed par t does not necessarily behave like the liquid-drop model 
function of deformation 



34j . See Ref. [35| for recent developments. Another method 



is to apply the semiclassical Wigner-Kirkwood expansion of the single-particle partition 

The relation to the Strutinsky shell correction method was 



function 
discussed 
developed 




and the treatment of realistic potentials with the spin-orbit term was 



This method was recommended in Ref. 
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26] to obtain the smoothed 



energy unambiguously. See Ref. 40) for recent developments. However, to achieve the 
same accuracy as the conventional Strutinsky shell correction, one has to include up to 
the third order terms in h 2 . The lowest term is nothing but the Thomas- Fermi energy. 
The calculation is rather complicated especially for the case without spherical symmetry. 
It should be also noticed that the semiclassical level density diverges at the threshold (or 
the barrier top in the case with the Coulomb potential), which has non-negligible effects 



for drip line nuclei [26]. In this paper, we stick to the conventional energy smoothing 
procedure and do not consider these other possibilities. 

The second difficulty in the shell correction approach with the continuum states in- 
cluded is the treatment of the pairing correlation, for which the simple BCS approximation 
is usually used with the "diagonal" (seniority) pairing force. The force strength is fixed 
according to the model space employed by the smoothed pairing gap method jsl, liol ]. 
Since the pairing model space is taken to be within about the major shell spacing above 
and below the Fermi level, the same problem as that of the smoothed energy arises for 
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unstable nuclei, where the Fermi level is so close to the threshold that the unbound states 
enter into the model space. This is a serious problem because finite occupation proba- 
bilities of unbound states lead to the formation of "neutron gas" surrounding nucleus. A 
complete solution of this problem requires the coordinate-space HFB method 4lj. The 
pairing energy is also affected by the continuum states in such an uncontrollable way that 
it increases infinitely as more number of states are considered. It is a major obstacle to 



the unambiguous prediction of the drip line 



161 ] . We extend the Kruppa's prescription to 



the treatment of the pairing correlation and try to solve this problem. 

The third difficulty in the the shell correction approach, which is not particularly 
related to the unbound states, is the inconsistency between the Fermi energy of the chosen 
single-particle potential and that of the macroscopic part j^J ; this kind of problems do not 
appear in the selfconsistent mean field approach, since the single-particle potential adjusts 
itself to give the correct Fermi energy. Though this problem is negligible in stable nuclei, 
it becomes severer near the particle threshold, which easily shifts the drip line by about 
te. particle numb e, in R e t .Q, P — ofaWoods -Saxon potential are adjusted in 
accordance with the bulk nuclear asymmetry of the droplet model; it is found that a fine 
tuning is necessary to obtain the coincidence of the Fermi energy of the adjusted potential 
with that of the macroscopic part. In this paper we solve this problem with an automatic 
adjustment of the potential depth in the Thomas- Fermi approximation. 

The main purpose of the present work is to solve the difficulties of the conventional 
microscopic-macroscopic approach. We propose remedies to all the three difficulties men- 
tioned above. Although our remedies are not perfect ones, we believe that a combined use 
of them gives much more reliable results for the shell correction calculations of unstable 
nuclei near the drip lines. This article is organized as follows. In Sec. Ull the present status 
of the shell correction method is reviewed and its difficulties are discussed in details. A 
new interpretation of the Strutinsky energy smoothing is also given there. In Sec. IHH 
the solutions to the difficulties are presented and the qualities of the improvements are 
examined in detail. Sec. IIVI is devoted to the conclusion. 
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II. THE PRESENT STATE OF THE SHELL CORRECTION METHOD 



A. The Woods-Saxon potential 

The Woods-Saxon potential is a finite-depth potential, which has a continuum spec- 
trum of unlocalized states as well as a discrete spectrum of localized states. Combined 
with a spin-orbit and proton's Coulomb potentials, it resembles very well the potentials 
for a nucleon in atomic nuclei, having a flat central part and a short-range tail. The 
expression we employ is given by 

.2 1 



Hws = ^ + VcE + Vso + i (1 ~ rs) Vco ' 



where the central part Vce and the spin-orbit part Vso are the standard ones [43] 



(1) 
(2) 

^ 2m re dC ) v ' % 

where m re d = ~inW with m being the nucleon mass, T3 the third component of the 
nucleon's isospin multiplied by two (1 for neutrons and —1 for protons), cr the Pauli 
matrix for the nucleon's spin (s = § cr), and the function Vws is defined by 



so 



Aso 

A-l, 



h 



Vce — Vws(f; Voce, ^ce, -Roce, ^ce, 0), 
2 



VVws(t*; VocE,«so,-Roso,aso,^)) • (cr x - V 



Vws(r; Vo, «, i? , a, /3) = -V 



1 ± K- 



N - Z 



A 



exp[D(r;R ,f3)/a] 



(4) 



Here, N, Z, and A are the neutron, proton, and mass numbers, respectively, while ± in 
front of k means + for proton and — for neutron. D(r; R , f3) denotes the (perpendicular) 
distance (with minus sign if r is inside the surface) between a given point r and the nuclear 
surface, so that D(r; R , f3 = 0) = r — R for spherical shape. The surface is specified by 
the radius Rq and the deformation parameters f3 = (f3\) as, 



R(9;Ro,P) = Roc v {(3) 



I + E^AOW 



(5) 



where the constant c v (f3) takes care of the conservation of the volume inside the surface 
against deformation (c v = 1 for (3 = 0). We consider only axially symmetric deformations 
and take into account the quadrupole (/3 2 ) and hexadecapole (^4) ones in this paper. The 
Coulomb potential Vco acts only on protons, and is created by electric charge (Z — l)e 
distributed uniformly inside the nuclear surface given by Eq. (jSJ) with R = Rqce- 
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The parameters -Roce (-Roso) and ctcE («so) are the radius and the surface diffuseness 
of the central (spin-orbit) potential. For N = Z nuclei, the depth of the central potential 
is Voce, while a dimensionless parameter Aso specifies the depth of the spin-orbit potential 
relative to the central potential. The quantities kqe and nso describe the nuclear isospin 
dependence of the two potentials. Note that the radii and diffusenesses of the central and 
spin-orbit potentials are different in general, but the shape of nuclear surfaces are taken 
to be the same, i.e., the common deformation parameters f3 are used in both of them. 
The set of the values of these parameters mainly used to obtain the results shown in this 
paper is the universal parameter set of Ref. 43J] (Note that the parameter ro-so(P) = 1-20 
in Table 1 of Ref. 43] is a misprint and should be replaced to 1.320, see Ref. }44|). It 
should, however, be noted that we modify the depth of the central potential in order to 
be consistent with the liquid-drop Fermi energy; see Sec. IIII HI for details. 

The Nilsson potential is a harmonic oscillator potential combined with a spin-orbit 
and an I 2 terms. Since its depth (or height) is infinite, its spectrum does not have a 
continuum part. See, e.g., Ref. [45J for the equations to define the potential. We employ 
the N osc - dependent Is and I 2 parameters of Ref. 46]. 

We use the anisotropic harmonic oscillator basis to diagonalize these single-particle 
Hamiltonians. The oscillator frequencies, U3 and u±, along the symmetry axis and 
the perpendicular axis, respectively, are determined by the two conditions; the volume 
conservation and the condition that they are inversely proportional to the root mean 
square length of each axis, which is calculated assuming the uniform sharp-cut density 
inside the nuclear surface given by Eq. ([51). Namely the conditions are u^uj_ = Uq and 
W3/w_l = \J (x 2 )uni/ {z 2 )uni , where ujq is the frequency for spherical shape and ( ) un i de- 
notes average value based on the uniform sharp-cut density. The number of the basis 
states is specified by the total oscillator quantum number N osc = n± + n 3 (n± = rii +n 2 ), 
where n« (i — 1,2, 3) is the number of oscillator quanta in the i-th axis. In the following 
discussions, we use the standard harmonic oscillator energy, hu = 41 /v4 1/3 MeV, and the 
Woods-Saxon potential is diagonalized in the oscillator basis with a frequency uq = 1.2 uj. 
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B. The shell correction method 



In the shell correction method, the total energy of a nucleus is assumed to be decom- 
posed as 

E = E mac + £ (E^ + E&) , (6) 

q=n,p 

where E mac is the energy of a macroscopic model like the liquid-drop model while E^ and 
£?pai r are the microscopic corrections. Because the equations to define the contributions 
from neutrons (q=n) and protons (q=p) are very similar, we show only the terms for 
neutrons in the rest of this paper. For the sake of conciseness, we omit the superscript 
(n) for the most part, i.e., E sh and -E pa i r designate E^ and E^ ivJ respectively. 
The term E^ is the shell correction energy, which is defined by 

E s h = E sp , — E'g.p.. (7) 

The first term on the right-hand side is the sum of the single-particle energies of occupied 
levels, 

N 

E s . p . = ^2u, (8) 

where e\ < 62 < • • • are the neutron single-particle energies. Since we are going to discuss 
about the Kruppa method (see Sec. IIIDI) . these levels are the results of diagonalizations 
of the single-particle Hamiltonian in a truncated harmonic oscillator basis and thus they 
are discrete through negative and positive energies. 
By introducing the (single-particle) level density, 

0(<O = (9) 

i 

the quantity E B , P , in Eq. (jSJ) can be written as an integral, 

E s . p . = f X eg(e)de, (10) 

J — 00 

up to the Fermi energy A. Analogously, the second term on the right-hand side of Eq. (j7]) 
is the integral of the product of the energy and a smoothed level density g(e) over a 
semiinfinite energy interval up to the corresponding Fermi energenergy A, 

E s . p . = I* eg(e)de, (11) 



with A determined to satisfy a constraint on the number of particles, 

[ X ~g(e)de = N. (12) 
The term E pa ,i T is the correction for the pairing energy, which is defined by 

-Epair = (-EbCS ~ #s.p.) — (-^BCS ~ ^s.p.) • (13) 

-Ebcs and .Ebcs are the energies of the BCS solutions of the pairing Hamiltonian with 
discrete and smoothed level densities, respectively. The terms in the first parentheses in 
the right-hand side represent the energy gain due to the pairing correlation. The terms 
in the second parentheses are the part of the pairing energy gain smoothly changing as 
a function of N and Z, which should be subtracted since it is already included in E mac . 
Explicit expressions for these quantities are given in Sees. IIIIDI and IIIIFI 
Using Eqs. (j7|) and ( 1T5|) . one can simplify Eq. (jBJ) as 

E = E mac + ( e bcs - e bcs) ■ ( 14 ) 

q=n,p 

However, from a physical point of view, we discuss E s ^ and E pa i T separately. It may be 



worth noticing that one often uses simplified expressions for t 



re smoothed part of the 



pairing energy in many of existing calculations, e.g., Refs. |l_0|, Il5| . assuming that the 
single-particle levels are uniformly distributed with the smoothed level density at the 
Fermi energy. In such cases Eq. (I14p does not hold exactly. In this paper we calculate 



-^bcs consistently without such simplifications, as is discussed in Sees. IIIIDI and IIHF 
As for the energy of the macroscopic part, we use the liquid-drop model of Ref. 47] 
this paper; see also Ref. |48] for its explicit form. 



in 



C. The Strutinsky smoothing method as a low-pass filter 

In the conventional Strutinsky smoothing method, the smoothed level density is ob- 
tained by a convolution integral with respect to the single-particle energy, 

m = - r 9W)f p (—) de' } (15) 



7 J -co \ 7 
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/CO 
f p (x)dx = 1, while 7 is the width 
-00 

parameter. The smoothing function is chosen as 

f p (x) = -^e^L^ix 2 ). (16) 
V 71 " 

Here, Ll/ 2 (x) is a polynomial of order p (the generalized or associated Laguerre polyno- 
mial [49]), with which the transformation ( jl5j) leaves g(e) unchanged, i.e., g(e) = g(e), if 
g(e) is a polynomial of order 2p. Note that the order of polynomial is denoted by "p" 



in, e.g., Refs. 



16 



26 



27| , so that the parameter p in these references is 2p in this work. 



Fig. [T] (a) shows the graphs of f p (x) for several values of p. 
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FIG. 1: The Strutinsky smoothing function in panel (a) and its Fourier transform in (b). The 
parameter p is a half of the order of the polynomial part of the function. 



For a discrete level density, 



9^) = I>(e-ei), 



the smoothed density is given by 



8=1 



M 



e - £i 



(17) 



;is) 
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where M is the number of single-particle levels included in the calculations. Owing to 
the gaussian form factor in f p (x), this transformation has a short-range character, which 
is a large advantage because it makes high positive energy levels unnecessary to evaluate 
Eq. (llip since they hardly affect g(e) at negative energies. 

Let us unveil another aspect of this transformation. The Fourier transform of a convo- 
lution of two functions is proportional to the product of each function's Fourier transform. 
Therefore, by denoting the Fourier transform of a function F as F like 

F{k) = / F(x) e~ ikx dx, (19) 

one can rewrite Eq. ffl5|) as 

9{r) = Uir)g{r). (20) 

(A similar expression in terms of the Laplace transformation is given in Ref. 38] in a 
different context.) The "wavenumber" r in Eq. ( 1201) has a dimension of (energy) -1 and 
may be regarded as a time variable (divided by h). Now, we show that the function f p 
has a typical shape of a low-pass filter. The Laguerre polynomial can be expressed in 
terms of the Hermite polynomials H 2 i as, 

^ V) = J2 QH 2l (x), Q = (~Y(2 2 H\)- 1 . (21) 

1=0 

2 

By multiplying e~ x to the generating function of Hermite polynomials, 

oo n 

-s 2 +2xs 



E#n(*)^> ( 22 ) 

one obtains 



n=0 U - 



e 



~{S-Xf _ V- TT fS -X 2 ' S 



J2H n (x)e-^°-. (23) 



n=0 U - 



The Fourier transform of the left-hand side can be calculated easily as 

/OO 00 q n 

e _ (s _ x) 2 e -ikx dx = ^ e - k */i e -iks = y 0F(_ifc) n e -* a /4* (24) 

which means that the Fourier transform of H n (x)e~ x2 is 

[°° H n (x)e~ x2 dx = v^F(-^) n e- fc2/4 . (25) 

J — oo 
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Using above results, one obtains the Fourier transform of the Strutinsky smoothing func- 
tion as 



p 1 (k\ 21 



e-^\ (26) 



1=0 

The term in the square brackets is the Taylor series of e^l 2 ^ 1 truncated at order 2p. For 
k <C kp Ut = 2y/p, the term is very close to e^ fc / 2 ^ and hence f p ~ 1, i.e., the filter is 
almost perfectly transparent. From this fact, one may give an alternative definition of the 
polynomial part of the Strutinsky smoothing function: It is a polynomial which minimizes 
the distortion of this low-pass filter in such a way that f^ l \0) = for 1 < I < 2p+ 1. For 
k 3> kp Ut , the term in the square brackets is much smaller than e^l 2 ^ 2 and thus f p ~ 0, 
i.e., the filter is almost completely opaque. 

TABLE I: Changes in the characteristics of the low-pass filter f p , i.e., the Fourier transform of 
the Strutinsky smoothing function, versus the order p of its polynomial part. The normalization 



1- fp 
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inverse funct 
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14.236 
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2.561 


100 


20.067 




2.007 
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In Fig.[T](b), the Fourier transform f p (k) of the smoothing function is shown for several 
values of p. One sees that they are almost a constant function near k = and decrease 
monotonically to zero. They become a half of the maximum around k w fc£ ut = 2^p 
(except p=0). The length of the interval where the function drops from 90% to 10% of 
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the maximum (/ p (0) = 1) is ~ 2.5 and almost independent of p. One can verify these 
features in Table [H In this way, the usage of higher order polynomials lengthens the 
width of the filter. At the same time, it shortens the width of the smoothing function in 
Fig. Q] (a) in a complementary manner. 
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^ 0.2 
0.0 
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1.0 
0.8 
^0.6 

0.2 
0.0 
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k' 

FIG. 2: Rescaled Strutinsky smoothing function in panel (a) and its Fourier transform in (b). 

Since the width of the filter of order p is proportional to v /p, it is convenient to use 
variables scaled with y/p, k! = k/*Jp. In Fig. [2] (b), we show f p {k'^/p) versus k! for 
several values of p. As the order p is increased, the cutoff becomes sharper while the 
position of the cutoff converges to k' = 2 independent of p; it approaches a step function 
9 (2 — \k'\) in the limit of p — > oo. Corresponding changes in the function f p can be found 
by using a dimensionless variable x' = v/p(e — e')/7 and a rescaled smoothing function 
fp( x '/y/p)/y/P to rewrite Eq. ( JIB]) as 

/oo 
g(e-'ya//y/p)fp(a//y/p)da//y/p. (27) 
-oo 

In Fig. [2] (a), f p (x'/ yjp ) / yjp is shown as a function of x' for several values of p. Although 
the convergence is slow, for very large values of p, fp{x'/^/p~)/^/p — (sin 2x')/irx', which 
can be obtained as the inverse Fourier transform of 9 (2 — \k'\). The curve for p = 50 is 
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quite close to this function in the interval shown in the figure. For larger x', however, the 
rescaled smoothing function decreases much faster than x'~ x due to the Gaussian form 
factor. 

Using the semiclassical periodic-orbit theory, the quantum mechanical level density ( jTTl) 
for a certain class of potentials can be represented by a sum of contributions from classical 
periodic (closed) orbits, the so-called trace formula |50l-l52l] . It is discussed that the origin 
of the gross shell structure can be well understood in terms of a few important short peri- 

*° rbitS1 f ° reXamPle ' beating pattejn of the level de^ in the spherical biUiart Q, 
the shell structure in deformed nuclei [54J , and the supershells in metal clusters [55( . The 
smooth part of the energy corresponds to the gross shell structure, to which only short 
periodic orbits contribute. The low-pass filter expression ( 120|) demonstrates clearly that 
the conventional Strutinsky smoothing cuts off the contributions of long periodic orbits 
with period (divided by ft) r = fc/7 3> Tf ut for the filter f p (k), where the cutoff period is 

r cut _ fccut^ = 2 ^ /7 _ 

If one changes 7 as 7 oc yfp for different choices of the order p in the smoothing 
function, the cutoff period r™* is independent of p, while the cutoff becomes sharper for 
larger p value as is clearly seen in Fig. |2] (b). In Sec. Ill Fl and Sec. IIII CI we use this 
fact for discussions on the plateau condition, i.e., the stability of the smoothed energy 
with respect to the smoothing width 7 and the order p specifying the curvature correction 
polynomials. 



D. Kruppa's prescription for the positive energy levels 

For finite-depth potentials like the Woods-Saxon one, positive energy levels appear as 
continuum states. They also affect the energy of bound nuclei through Eqs. ( JTTT) and (TI8I) . 
Their contribution becomes larger when A is closer to zero. If one obtains the positive 
energy spectrum by diagonalizing the Hamiltonian in a truncated oscillator basis, the 
positive energy spectrum is not continuous but discrete. Thus one can calculate the 
summation in Eq. ( JT8j) straightforwardly. However, the result depends strongly on the 
size of the basis M. In fact, the smoothed level density ([TBI) diverges in the continuous 
limit and so does the smoothed energy ( fill) ; it is monotonically increasing or decreasing 
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as increasing the size of the basis [if], • A practica 



way to avoid this problem is to 



15 



13| to take N™ x « 12 for the 



restrict the size of the basis; it is recommended in Ref. 
harmonic oscillator basis. However, in such small bases, negative energy levels may not 
be sufficiently accurate as we will see in the followings (see, e.g., Fig. HJ). 

A way to circumvent the diverging single-particle level density g(e) due to the particle 
continuum is to replace it with the so-called continuum level density g c (e). In the case of 



spherically symmetric potentials, it is written as 
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20|, 



i:bound Ij 



(28) 



where 5y(e) is the scattering phase shift. This expression was used for the calculation of 
shell correction energy and it was found that the contributions of the particle continuum 
(the second term on the right-hand side) through E s . p . are never negligible even in stable 



nuclei for finite-depth potentials 18j, |21|. 



One can roughly regard the continuum level density as the difference between the full 
and the free level densities 20] . Taking the energy derivative of the phase shift in Eq. (1281) 
means calculating the level density from the number of states. The number of states is 
actually proportional to the phase of the radial oscillation of the wavefunction because 
an increase in the phase by tt corresponds to the addition of one radial node in the box 
boundary condition. The phase shift is the difference of the phases between full and free 
solutions. Therefore, the definition in terms of the phase shifts is actually equal to taking 
the limit of infinite volume of the difference between the full and free level densities in a 
finite volume cavity. 



This can be shown more rigorous 



symmetric potentials is given by 
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y. The generalization of Eq. ( 128]) for non-spherically 
23] 



1 



tr f 



dS(e) 

'-3T 



(29) 



where S(e) is the on-shell S-matrix corresponding to the single-particle Hamiltonian H 
with energy e, and tr e means the restricted trace operation with respect to the eigenstates 
with energy e. Note that Eq. f]2"9"]) contains the contributions from the bound states because 
they appear as poles of the S-matrix. This quantity is related to the time-delay [56( , and 
shown to be identical to the trace of the difference between the full and free single-particle 
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Green's functions 23J]. In this way, the level density can be written as 

(30) 



g c (e) = -Im 

7T 



1 1 

tr tr ■ 



H-e H -e. 

where H is the free Hamiltonian (with the repulsive Coulomb potential for proton), and 
tr here is the full trace operation. This expression clearly tells that both full and free 
level densities are divergent for positive energies but their difference is finite. It is used for 



investigation of the level density in Ref. 24| by using the Green's function technique 



3- 



Inspired by Eq. f[3"Uj) . Kruppa has introduced a prescription 22J, which is suitable 
to treat the particle continuum by the diagonalization method with, e.g., the harmonic 
oscillator basis. He has demonstrated that results with his prescription have much weaker 
dependence on the size of the basis and converge for enough large basis. Let us call 
his prescription the Kruppa method. This method changes the definition of g(e) as the 
difference of the single-particle level density between the full Hamiltonian (including the 
potential) and the free-particle Hamiltonian, 

M M 

9(e) =► g K (e) = £<5(e-e,) -£<5(e-e°), (31) 

i=l j=l 

where and e° are the eigenvalues of the full and the free Hamiltonians, respectively. 
Here M is the dimension of the basis commonly used in the two diagonalizations, and 
g K (e) — >■ g c {e) as M — > oo (see Eq. (1301 ). Note that, for one-body observables like the total 
single-particle energy in Eq. (flUj) . the free energy terms in Eq. (|3T]) do not contribute as 
long as A < 0, i.e., when the Fermi energy does not exceed the particle threshold. However, 
they contribute to the smoothed quantities. Now the smoothed level density g(e) should 
be obtained by applying the Strutinsky smoothing to this g K (e): 

9(e) ~9 K (e) = -j:fJ^)--j:f P ( e -^-). (32) 

7^ V 7 J Ij^x V 7 J 

The redefined level density g K (e) converges to g c (e) for sufficiently large basis sizes, the 
reason of which is explained transparently in Sec. Ill El 

The continuum level density was originally used to calculate the second virial coefficient 
(related to the deviation of the equation of state from that for the ideal gas) arising from 

n □ 

the interaction between gas particles [13, [20] . For this purpose, one naturally has to 
separate the part corresponding to the free motion of non-interacting particles from the 

16 



integral over the continuous spectrum. Unlike this case, the reason to subtract the free 
spectrum is not so obvious in the calculation of the shell correction. At present, we do 
not know whether it can be derived rigorously from a more basic theoretical framework. 
Nevertheless, it certainly seems to be the most reasonable prescription so far to obtain 
physically meaningful results. 
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FIG. 3: Smoothed level densities for the full and free Hamiltonians and their difference (i.e., 
the Kruppa's level density) as functions of the single-particle energy in MeV. The smoothing 
parameters used are 7 = 1.2 hu and p = 3. Panels (a), (b) and (c) are for N™ x = 12, 20 and 
30, respectively. The nucleus is 164 Er with deformation fii = 0.27 and $4 = 0.02. 



In Fig. [3] we show three kinds of level densities, i.e., the full (with Woods-Saxon poten- 
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tial), the free, and the Kruppa for 164 Er. They are the results of the Strutinsky smoothing 
with 7 = 1.2 hu) and p = 3. The potential is deformed with /3 2 = 0.27 and /?4 = 0.02. 
The number of basis states M in Eq. (13T]) is specified by the maximum number of the 
oscillator quanta A"™ c ax , M = \(N£%* + l)(N™ x + 2)(N™ X + 3). Comparing panels (a), 
(b) and (c), one can see that positive energy part of the full and the free level densities 
are increased rapidly as N™ x is increased from 12 to 20 and to 30, while the Kruppa's 
level density does not change essentially. This clearly shows the fact that continuum parts 
of both the full and free densities are divergent but their difference is convergent. The 
energy range of the most influential part is e < A for the smoothed single-particle energy 
E B .p. and |e — A| < A ~ hu> for the smoothed BCS energy E^cs- Though the difference in 
this part between the calculated level density with N™ x = 12, 20 and 30 is much smaller 
than that in positive energy, e.g., at e ~ lOMeV, it brings about significant differences to 
the resulting nuclear properties, especially to the pairing correlation (see Sec. IIII Gj) . 

All the smoothed quantities in the Kruppa method are obtained by replacing q(e) with 
<? K (e). The shell correction energy E s ^ by this prescription is investigated in Ref. [27}, and 
shown to be also convergent when increasing the size of the basis. Examples are depicted 
in Fig. H] as functions of the basis cutoff N™ x - Without the Kruppa's prescription, 
depends on the size of the basis even in a stable nuclei 166 Er, and the dependence is much 
stronger in a neutron rich nuclei 226 Er. The subtraction of the continuum contributions 
reduces the dependence on the model space drastically and the shell correction energy 
with the Kruppa method converges in the large N™ x limit. These examples clearly show 
that the Kruppa's prescription is indeed promising. We extend it for other observables in 
SecIIlEl 

It is also worth noting that while the shell correction energy E s ^ almost converges at 
N™ x ~ 12, the sum of the single-particle energies E s p . itself does not; especially for the 
unstable nuclei 226 Er the single-particle energies are not obtained very accurately when 
A^™c X — 20. From this fact one may compose a syllogism on the necessity of the Kruppa 
method. (1) For large N™ x , the Kruppa method is necessary to treat correctly the dense 
positive energy spectrum. For small N™ x , it is not necessary. (2) One has to use large 
N™ x for sufficiently accurate bound-state energies. (3) One needs the Kruppa method. 
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FIG. 4: Neutron shell correction energies E s h ((a),(c)) and the sum of single-particle energies 
-^s.p. ((b), (d)) for 166 Er ((a),(b)) and 226 Er ((c), (d)) as functions of the oscillator basis cutoff 
iV™c X . The values of E s ^ obtained with and without the Kruppa's prescription are designated 
by filled circles and squares, respectively. The deformation parameters are = 0.280 (0.255) 
and /3 4 = 0.005 (-0.033) for 166 Er ( 226 Er) while the smoothing parameters are 7 = 1.2 froo and 
p = 3. 

E. Oscillator-basis Thomas-Fermi approximation for Kruppa's level density 



One can roughly reproduce the shape of the Kruppa's level density in terms of a new 
variant of the Thomas-Fermi approximation within the limited phase space corresponding 
to the truncated oscillator basis. We call it the oscillator-basis Thomas-Fermi (OBTF) 
approximation in this paper. One can also demonstrate the independence of the results 
of the Kruppa method from iV™^ x (if it is sufficiently large) in this approximation. 

We study only spherically symmetric potentials without spin-orbit couplings. Lifting 
these restrictions is possible but does not seem to be very meaningful, because it turns 
out that the OBTF approximation is not sufficiently accurate to be used as a replacement 
of the smoothed energy in the realistic Strutinsky calculations. This corresponds to the 
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known fact that, in order to obtain the Strutinsky smoothed energy, one has to include up 



38|, 



, in which the Thomas- 



to the third order terms in the semiclassical h 2 expansion 
Fermi approximation is the lowest. 

Hence we express the Hamiltonian for a nucleon as 

H (P' r ) = ^ + V ( r )i Vr ( r ) = VA c E (r) + ^(l-r 3 )K ; o(r), (33) 

where VcE( r ) and Vco( r ) are the central and Coulomb potentials in Sec JII Al with spherical 
shape, i.e., with the deformation parameters f3 = 0. The states are assumed to be doubly 
degenerated for the two spin states s z = ±|. In the Thomas- Fermi approximation, the 
number of particles in the potential well for a given single-particle energy e is given by 



o 



oo 

4tt / p TF (r,e)r 2 dr (34) 



where p TF (r, e) is the particle density at position r for Fermi level e expressed as (using 
the Heaviside function 9), 

(2m) 3 / 2 



Ptf\T, e J 



37T 2 ft 3 

The level density is related to the number of particles T(e) as 



-V(r)\ 3/2 6(e-V(r)). (35) 



dT(e) roc dp (r,e) r2dr ^ 



TF de Jo de 

This level density diverges above the particle threshold, e > 0, for finite-depth potentials 
because of the infinite volume of the space. 

The idea of OBTF is to extend the Thomas- Fermi approximation in such a way that the 
phase space is limited within a subspace spanned by a truncated harmonic oscillator (HO) 
basis, which can be specified by the maximum kinetic energy as a function of position as 
in the followings. By replacing V(r) with the oscillator potential Vno( r ) — \muj 2 r 2 in 
Eq. (1351) . one obtains 

WeH^T^) 3 , (37) 



3 Xhuj , 

which is always finite. A truncated oscillator basis is usually defined by the maximum 
oscillator quantum number N™ x - Equating the right-hand sides of Eqs. ( 1371) to the 
number of states with A^ osc < N™ x leads to the cutoff energy of the truncated basis, 
e = e cut (N™*), 

e cut = hu [(AC X + l)(iVr x + 2)(iV™ x + 3)] 1/3 . (38) 
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We also define -R max by a condition Vno(Rn 

R 

max 



-cut ) 



i.e., 



' 2e cut 
mco 2 ' 



(39) 



and the local maximum kinetic energy expressed as 



-kin 



(r) = (ecut-%o(r))fl(e cut -F H o(r)) 



-moj 
2 



(R 



2 

max 



r 2 9 (R n 



Now we define the level density in OBTF, similarly to Eq. ( 136|) . as 



9c 



An 



(40) 



(41) 



io de 

where p OB (r, e) is equal to the right-hand side of Eq. (135]) with an additional restriction 
that the energy e should be smaller than e^ x (r) + V(r). Its derivative is given by (with 
the (5-function contribution from the Heaviside function having no effects), 

dp °f e ' £) = \e ~ V{r)\ 1 ' 2 9 (e - V{r)) 6 (e^(r) + V(r) - e) . (42) 

In this way, the finite level density <? OB (e) is obtained for a given maximum oscillator 
quantum number 



jmax 
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FIG. 5: A schematic figure to explain the convergence of the Kruppa's level density in the 
oscillator-basis Thomas-Fermi approximation. The abscissa represents the radius r from the 
center of the nucleus, while the ordinate is the kinetic energy ek m = e — V(r) of a single nucleon. 
Hatched area A, B, and C are the domain of integrations to obtain r OB , Fq B , and Tq B , respec- 
tively. Parabolas drawn with dash and dot lines stand for the maximum kinetic energy e{5n X ( e ) 
in the oscillator basis with N™^* and N™ xl (> N™ x ), respectively. 



The Kruppa's level density g^ B (e) is defined as g t 



the free-particle level density expressed as g® B (e) 



OB v 

4tt 



■5o B (e), where fi_[c) js 



^OB 



dp Q (r, e) / de r 2 dr with p° (r, e) 
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obtained by omitting V(r) in the right-hand side of Eq. (142]) . (This is for neutrons and 
changes necessary for protons are described in the following paragraph). It is readily 
shown that 

F f \ [" t <\A> 2 ( 2m ) 3/2 f 1/2 2, , ( A o\ 

r o B (e) = / 9 A e ) de = is — / e kin r de kin dr, (43) 

where ekin = e' — V(r), and the domain A of integration is depicted in panel (a) of Fig. El 
Changing the domain to B and C shown in panels (b) and (c) gives the similar expressions 
for T° B (e) = f g° B (e')de' and r£ B (e) = f g^ B (e')de', respectively. By enlarging the 

J -co J— oo 

harmonic oscillator basis (i.e., by increasing iV™^ x to N™ x> > ^^c X )' domains A and 
B are expanded while domain C is left unchanged. The unchanged domain results in an 
unchanged number of levels and thus an unchanged level density. This explains pictorially 
why the Kruppa's level density converges for large N™ x above the particle threshold, 
e > 0. 

It is also possible to show that ^ B (e) oc e -1 / 2 as e — > oo after the limit N™ x — > oo 
is taken. For an arbitrarily given e > 0, one can take sufficiently large N™ x to express 
the region C as {(e^n, r) | < r < oo, e < e^n < e — ^( r )} with an approximation that 
V(r) = for r > R max to obtain 

\3/2 



T K (e 



^r{le-V(r)r-^}^r. (44) 
67in Jo 1 > 



Thus, for the level density S'obI 6 ) = dT^ B (e)/de, 

1/2 k / \ 2(2m) 3/2 r°° V(r)r 2 dr 

6 ^ l + [l - V(r)/e]^ 

-> - V ' / V(r)r 2 dr (e ->• oo). (45) 
7m Jo 

It can be confirmed that the following expression is a very good approximation for the 
iVosc* — » oo limit of the Kruppa level density in the whole range of single-particle energy: 

^b( £ ) « 1 p / (e " ^(r)) 1/2 5 (e - 7(r)) - e 1/2 # (e) r 2 rfr. (46) 

For protons, one can repeat the same argument if one includes Vco( r ) i n the free 
Hamiltonian, because Vco( r ) is not negligible even at r = R max and V(r) includes Vco( r )- 
In the same way, it is readily seen that the Kruppa level density is convergent as N™ x 
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oo and in a very good approximation, 
2 (2m) 3 / 2 



'OB V 



( c _ V (r)f 2 9 (e - V(r)) - (e - V co (r)) 1/2 (9 (e - V co (r)) 
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FIG. 6: Neutron's level densities for the full and free Hamiltonians and their differences obtained 
in the oscillator-basis Thomas- Fermi approximation ((a),(c)) or with the Strutinsky smoothing 
method ((b), (d)). The smoothing parameters are 7 = 1.8hu and p = 3. The oscillator basis has 
N™ x = 12 ((a),(b)) or iV ™ c ax = 20 ((c), (d)). The nucleus is 154 Er. The potential is spherical 
(02 = Pa = 0) and the spin-orbit potential is turned off. 



In Figs. [6] and [7J the OBTF level density is compared with the smoothed exact level 
density for a spherical nucleus 154 Er. Figure [7] includes the proton Kruppa densities. The 
spin-orbit force is neglected and a larger smoothing parameter 7 = 1.8hoj is used with 
p = 3 for this calculation to make the comparison more appropriate. One can see that 
the OBTF is a fairly good approximation for both N™ x = 12 and N££* = 20. The 
proton Kruppa level density is very similar to the neutron one except that the single 
particle energy is shifted by the Coulomb barrier, about 10 MeV in this nucleus, as is 
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FIG. 7: Kruppa Level densities obtained in the OBTF approximation and with the Strutinsky 
method for 154 Er; panel (a) is for neutron and (b) for proton. The arrows denote the particle 
threshold (e = 0) . The same calculation as in Fig. [U] is used except that the oscillator basis has 



30. 



shown in Fig. [7J Threshold behaviors of the OBTF neutron and proton level densities 
are slightly different, which reflects the effect of the long-range Coulomb potential, while 
such differences do not exist for the smoothed exact level densities. 

Apart from oscillations at large e, the average behavior of the continuum level density 
is well reproduced in the OBTF approximation in Fig. [HI It can also be clearly seen how 
subtracting the free level density, Eq. (131 D . works to diminish the dependence on N™ x in 
the Kruppa method. However, as it is clearly shown in Fig. [7J the precise shape of the 
smoothed level density cannot be obtained; especially, the peak near the threshold (e ~ 0) 
shows cusp behavior in the OBTF density, which is characteristic to the semiclassical 
approximation, while the smoothed density looks like a broad peak. In order to obtain 
more precise level density one has to go beyond the Thomas-Fermi approximation. 

The Kruppa OBTF level density g~ (e) becomes negative in the range of single-particle 
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max 4-1, 
osc 5 Wlcll 



energy, e cut + V(0) < e < e cut (see Fig. [6] (a)), and it can be shown, for a finite N ( 

/oo 
g^ B (e)de = 0. This behav- 
-oo 



ior is also known 



24| in the exact continuum level density g c (e) defined by Eqs. (128]) — (1301 . 



/oo 
g c (e)de = 0. In this way the OBTF 
-oo 

Kruppa level density satisfies the desired property of the continuum level density. 



F. Plateau condition 

It would be preferable if E s _ p _ of Eq. ( ITT]) did not depend on the parameters concerning 
the smoothing of the level density (7 and p in Eq. (JT5i) ) because their values can be chosen 
arbitrarily. Since a perfect independence is unlikely to be satisfied, one usually demands 
a weaker condition that the dependence is very weak in a certain interval of 7 for a few 
values of p. This is the meaning of the plateau condition in this paper. 

For the Nilsson spectrum, a long plateau appears in most cases (see, e.g., Ref. 0, 



21] ). On the other hand, for finite-depth potentials like the Woods-Saxon potential, the 



situation is subtle. If the oscillator basis is truncated at iV™!™ m 10 to 12, reasonable 



plateau are obtained in many eases and .\'™' x = 12 is a. recommend value <is a. 



osc 



working prescription in Ref. [15j. However, the model space defined by N™ x — 12 
is not large enough to calculate single-particle states accurately, especially for unstable 
nuclei (see, e.g., Fig. UJ (d)), and this truncation is not justified. The appearance of 
plateau obtained by the relatively small model space with iV™^ x m 10 to 12 is accidental 
and increasing N™ x drastically change the situation [l^] ; the shell correction energy 
depends strongly on the smoothing width 7. This clearly indicates that a naive inclusion 
of continuum states by the diagonalization method does not work. Then, the continuum 



level density, Eq. (|28]1 . is used to calculate the shell correction energy jl8|, |2lJ. Although 
the dependence on 7 is weaker if the phase shift is calculated up to enough high energies, 
no good plateau like in the case of the Nilsson potential is obtained 2s| . 

We show examples in Figs. |8]to[TTJ Figures [8] and [9] depict the neutron shell correction 
energies E s ^ calculated with the standard Strutinsky smoothing method and with the 
Kruppa-method, respectively, changing the basis size specified by the maximum oscillator 
quantum number N™ x - The results for the stable and unstable nuclei, 166 Er and 226 Er, 
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FIG. 8: Neutron shell correction energies without using the Kruppa's prescription as functions 
of the smoothing parameter 7 in units of hco. Each curve represents the result with different 
order p = 3 to 6 of the smoothing function (|16p . The diagonalization basis is N™^* = 12 
((a),(b)), = 20 ((c),(d)), and iV™ x = 30 ((e), (f)). The nucleus is 166 Er ((a),(c),(e)) and 

226 Er ((b),(d),(f)). The deformation parameters are the same as in Fig. SI 



are compared. If N™ x — 12 is used for the basis size, a plateau-like behavior is observed 
in a reasonably long range for 166 Er and in a shorter range for 226 Er. But this is "spurious" 
because, using larger N^ x , the shell correction energy depends more strongly both on 
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FIG. 9: Same as Fig.[8]but with the Kruppa's prescription employed. The results with N ( 
20 are omitted since they are very similar to the results with N™ x = 30. 



max 
osc 



the smoothing width 7 and the order p of the curvature correction polynomial, while the 
range of "real" plateau in the case of the harmonic oscillator potential grows as the basis 
size increases 25[. The possible reason of this "spurious" plateau is that the number 
of discretized continuum states with iV™^ x = 12 is just suitable for the smoothed level 
density to be approximated by the lower order polynomial functions across the particle 
threshold e ~ 0. Increasing the basis size the curvature of the smoothed level density 
changes suddenly at e « 0, as is shown in Fig. [3j which no longer be approximated by a 
simple polynomial; leading to the strong dependence of E s ^ on 7 and p. Therefore, it is 
difficult to obtain reliable shell correction energies in the standard Strutinsky method. 

In contrast, the Kruppa's prescription reduces the basis-size dependence dramatically, 
as can be seen in Fig. [91 Compared with the standard method, where the plateau condition 
is more and more unsatisfied as increasing the basis size, the stability against the increase 
of N™ x is a very important feature of the Kruppa method. However, although there are 
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almost degenerate local minima with different order p's, the plateau is not well established 
generally. The situation is worse for unstable nucleus 226 Er. A possible improvement will 
be discussed in Sees. II II Al to II II CI 




0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 




At first sight, the dependence of the shell correction energy on the smoothing width 
7 is quite different when the order p of the smoothing function (I16p is changed. Close 
inspection reveals, however, that the different curves in each panel of Figs. |S] and M are 
almost isomorphic if they are drawn as functions of the y/p -scaled width parameter, 

1p = 7/V5A (48) 

as shown in Figs. [10] and [TTJ Here we choose j p=3 = 7 because p = 3 is a standard choice 
for the curvature correction polynomial. The reason of this "isomorphism" between the 
results with different order p's can be understood from the discussion in Sec lII CI the 
range of the low-pass filter increases when employing the larger order p, and if is used the 
variable scaled with ^Jp the cutoff ranges are the same but the filter becomes sharper, 

28 




FIG. 11: Same as Fig. [9] but plotted as functions of the scaled smoothing width parameter 
7 P = 7/ y/p/3. Only the results of N™ x = 12 are shown since the results with different N™ x 
look similar. (See the results with = 30 in panels (a) and (c) of Fig. [191 ) 



as is clearly shown in Fig. [2] (b). Therefore, the complete isomorphism means the shell 
correction energy is independent of the sharpness of the filter. We have found that, for 
calculation with larger N™ x , better isomorphism is generally obtained by the Kruppa 
smoothing method than by the standard one; compare Fig. [10] with Figs. [HI Even better 
isomorphism is obtained in the improved treatment in Sec. IIII Al (see Fig. fT9l) . In the 
following discussions for the plateau condition, we always use the Kruppa prescription 
and show the results as functions of the y/p-scaled width parameter 7 p (j48j) . 

In the course of writing the present paper, we noticed that a similar scaled smoothing 



width is used for investigating the plateau condition in Ref. [58 



where the isomorphism 



of the smoothing width dependence between different order p's is not as good as in our 
case. This is due to a different choice of basic smoothing function that is not gaussian. 
We believe that the Fourier transform of the smoothing function will be useful for more 
detailed comparison of our results with those of Ref. [58(. 



G. The reason for no good plateaux 

A clue to find the origin of this difference between the Nilsson (or the harmonic oscilla- 
tor) potential and the Woods-Saxon (or the finite-depth, in general) potential is the fact 
that the Strutinsky smoothing is a low-pass filter as discussed in Sec. Ill CI In the Fourier 
transformed world, the smoothed level density is simply the original density multiplied 
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by the filter. Therefore, the Fourier transform of the original level density (fT7|) . or the 
Kruppa density fl3TT) . 



M 



9\r 



M M 

) = ^exp(-zrei) - ^exp(-ire 



(49) 



i=i i=i 

should be investigated. 

In this subsection, we employ the units hu for the energy (e, 7, and a) and (fro;) -1 for 
the Fourier transformed time variable r, and regard e and r as if they were dimensionless. 
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e r 

FIG. 12: Neutron's level density of the Nilsson Hamiltonian in panel (a) and the absolute value 
of its Fourier transform in (b) for 208 Pb. N™ x = 9 is used for the basis size. In panel (a), 
smoothed level densities with 7 = 1.2 fiu and with 7 = 0.5 fioj (p = 3 for both) are included, 
while the Fourier transform in (b) has not been smoothed and calculated directly by Eq. (I49p . 
The units are huj for e, (hu)^ 1 for r and g, and g is dimensionless. 



The level density g(e) and its Fourier transform g{r) for the Nilsson potential are shown 
in Fig. [12] for the spherical nucleus 208 Pb. The single-particle states up to N™ x — 9 are 
included because the Is and I 2 parameters are given only for them 46j] . Fig. [12] (a) shows 
that the smoothed level density with 7 = 1.2 Tluj is approximately a quadratic function in 
e, while the major shell oscillation is clearly seen in that with 7 = 0.5 fiu. This indicates 
that the semiclassical property of the Nilsson spectra is essentially the same as that of the 



HO potential; its Thomas-Fermi level density is g, 



HO/ 



e 2 / \huf, see Eq. (fSTJ). Note 



that the so-called "iso-stretching" is done for the neutron and proton frequencies in the 
Nilsson potential j45|, cu n = (2N/A) 1 / 3 u and u p = (2Z/A) 1 / 3 u, so that the coefficient of 
e 2 is reduced by a factor 208/(2 x 126) in Fig. QJ (a). 

Concerning the behavior of g(r) shown in Fig. [12] (b), one can see a very low density 
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interval (2 < t < 5). One may probably call its origin as the harmonicity of the potential. 
If the cutoff period r™* = 2 v /p7 _1 of the filter f p (rj) is in this interval, the result of the 
filtering, g(r), hardly depends on rf ut (see Sec. Ill Cj) . Namely, the dependences on 7 and 
p are weak and there appears a plateau. This feature can be qualitatively understood by 
considering the case of the anisotropic harmonic oscillator potential, for which the spectra 
are equidistant and the sum in Eq. ( 1^91) with M — > 00 can be evaluated as an infinite 
geometric series to be 

-1 



(2i) 3 sin ^—rhux^ sin (^rhuo^j sin (^—rhu : 



(50) 



with bj x w ojy ~ uj z ~ w, which has a long low density interval between r = and r = 27r 
(in units of {hu)^ 1 ). 

In contrast to the case of Nilsson potential, the Fourier transform of the Woods-Saxon 
level density shown in Fig. [13] does not have this low-amplitude region, although the 
smoothed level density in the panels (a) and (b) clearly shows the similar major shell 
oscillation to that in Fig. [T^] (a). In Fig. [12] the Fourier transforms of not only the Woods- 
Saxon spectra but of the Kruppa spectra and of the restricted spectra within the bound 
states are also depicted. The Fourier transform of the bound-states spectra is smaller than 
that of the Woods-Saxon spectra, but is still about a factor two to three larger than that 
of the Nilsson spectra in the low density interval (2 < r < 5) (note the difference of scale 
in ordinates in Fig. [121 and [TBI . Moreover, the Fourier transform of the Kruppa spectra 
is larger than that of the Woods-Saxon spectra on average, and increasing the basis size 
makes the situation worse. This clearly shows that the Kruppa prescription does not 
help to make a plateau in the shell correction energy, which is already confirmed in the 
previous subsection. The mechanism to develop a long plateau in the Nilsson spectrum 
is not functioning in the Woods-Saxon spectrum. This seems to be the very reason for 
the absence of plateau for the Woods-Saxon spectrum. 

The Woods-Saxon potential is different from the Nilsson potential not only in the 
anharmonicity but also in the finite depth. It seems interesting to investigate further the 
difference between negative and positive parts of the spectrum. In order to examine this 
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FIG. 13: Neutron's level densities of the Woods-Saxon Hamiltonian ((a),(b)) and their Fourier 
transform ((c), (d)). The curves in panels (a) and (b) are the results of smoothing with 7 = 0.5 Tiu 
and p = 3, while those in (c) and (d) have not been smoothed. The harmonic oscillator basis 
for the diagonalization is N™ x = 12 in panel (a) and (c), and iV™^ x = 20 in (b) and (d). The 
full and Kruppa spectra are shown with solid and dash curves, respectively. In panels (c) and 
(d), the Fourier transform of only negative energy levels are also shown with dot curves. The 
units are the same as in Fig. [12j 



point, the short-time Fourier transform 59] seems useful. It is denned by 



F(k, x\ a) 



2tt 



F(x')w a (x' - x)e~ ikx 'dx\ 

F{k')w a {k - k')e lk ' x dk', 



-ikx 



(51) 



where for the window function and its Fourier transform we employ 



W a (K) 



na e 



(52) 



In this paper we apply it to a "short energy interval" Fourier transform. Figures [T^l to [T^l 
show g(r, e; a) with a = \/2, which gives the same size of window widths in two comple- 
mentary variables r and e. The same nucleus 208 Pb is used for this calculation. The input 
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level density g(e) has been smoothed with 7 = 1.2 and p = 3 in the right-hand panels (no 
smoothing has been done for the left-hand panels). The location of cutoff due to these 
smoothings are r™* = 2 y Jp r )~ l = 2.9 for 7 = 1.2, but then the cutoff result is blurred by 
the convolution with the window function of width y2/a = 1 (see Eq. ( 15TI) ). 
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(d) Nilsson, N<=9, 7=1.2 
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FIG. 14: Absolute value of the short-time Fourier transform of the level densities of the harmonic 
oscillator ((a),(b)) and the Nilsson ((c), (d)) Hamiltonians. The units are hu for e and (huj)^ 1 
for r. See text for explanations. 



As for the harmonic oscillator spectrum shown in Fig. (a), there persist two main 
components r = and r = 2n irrespectively of e. The component r = 2n corresponds 
to the major shell spacing (hu). There is a large low-amplitude region between the two 
hills along lines r = and r = 2tc. The cutoff for the standard smoothing parameters 
7 = 1.2 and p = 3 is rf ut = 2.9, which is almost at the center of this region. The result of 
this standard smoothing is shown in Fig. (b), in which the hill at r ~ 2tt is removed 
completely while that at r ~ is left almost intact. This explains the existence of a 
perfect plateau. 

In the case of the Nilsson spectrum shown in Fig. [TD (c), the hill at r ~ 2tt becomes 
distorted and lowered, which reflects the disturbance that the spin-orbit and the I 2 terms 



33 



(a) WS, N<=12, 7=0 



(d) WS,N<=12, y=1.2 




6 
4 
2 

-2 
-4 
-6 



6 
4 
2 

-2 
-4 
-6 



2 4 6 8 
(b) free, N<=12, y=0 
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(e) free, N<=12, 7=1.2 
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(c) Kruppa, N<=12, 7=0 
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(f) Kruppa, N<= 12, y=1.2 
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FIG. 15: Absolute value of the short-time Fourier transform of the level density of the Woods- 
Saxon Hamiltonian calculated with N™ x = 12. The units are Tiio for e and (ftu;) -1 for r. See 
text for explanations. 

of the Nilsson potential bring to the periodicity with the major shell spacing. However, 
the hill at r ~ and the low-amplitude region between the two hills are almost the same 
as in the harmonic oscillator case. This clearly explains that the similar good plateau can 
be expected in the Nilsson potential. 

For the Woods-Saxon spectrum, we show the absolute values of the short-time Fourier 
transforms of the full ((a),(d)) and the free ((b),(e)) spectra as well as the absolute 
value of their difference (the Kruppa's level density) ((c), (f)) in Fig. [15] for N™ x — 12 
and in Fig. [16] for N™ x — 20. One can still see the valley between the two hills in 
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(b) free, N<=20, y=0 
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(c) Kruppa, N<=20, y=0 
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FIG. 16: The same as in Fig. [15] but calculated with N ( 
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the full (a) and the Kruppa (c) spectra at negative e. However, it is much more filled 
than for the Nilsson spectrum. At positive e, the landscape is too complicated to be 
regarded as a single valley. In the results of the standard smoothing in panels (d), (e), 
and (f), the contours are much more irregular than those in harmonic oscillator and Nilsson 
cases. These irregularities indicate the existence of nonvanishing structures grown in the 
valley. Because their contributions change sensitively by small shifts in the cutoff from 
the standard value, the plateau can be destroyed completely. 

In the meantime, comparing the Kruppa spectrum with iV™^ x = 12 and N™ x — 20, 
one sees that the spectrum does not change at negative energies but continues to change 
at positive energies versus N™ x - The changes at positive energies originate in both the 
full and the free spectrum. These time structures at positive energies are most likely to 
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be the remnant of the property of the diagonalization basis. 

A related fact is that there are no clear changes of the principal time component of 
the full Hamiltonian between positive and negative values of the single-particle energy e. 
One can see only obscure and N™* dependent changes. Such clear changes would occur 
if the major shell interval were changed altogether at e = 0. Indeed, in Appendix C of 
Ref. |6oj], Magner et al. seem to have obtained such a clear change in the major shell 
interval between negative and positive energies that they could remove the difference 
through a transformation of the energy to obtain a plateau behavior. The difference of 
the results between them and us seems to be originated mainly in the difference between 
solutions in an infinite wall and those in an oscillator-basis expansion. In the latter case, 
the shell structure at positive energies is thought to be strongly connected with the basis. 

III. IMPROVEMENTS TO THE SHELL CORRECTION METHOD 
A. Reference density method 

Although the dependence of the results on the smoothing width is unremovable com- 
pletely, it is still preferable to make it as small as possible. In the Kruppa method, this 
dependence comes principally from the diffusion of the peak of the level density at thresh- 
old energy (e ~ 0). This peak is so sharp that it is inevitably more diffused by larger 
widths. Since this peak exists already in the (oscillator-basis) Thomas-Fermi approxima- 
tion, it should not be diffused but be kept unchanged. 

We now propose a prescription to prevent this diffusion, which we call the reference den- 
sity (Strutinsky) method. In the method, one applies the Strutinsky smoothing procedure 
not directly to the original discrete spectrum but to its deviation from some continuous 
reference level density g re f, i-e., to g(e) — g rc f(e). Note that our concern is the Kruppa level 
density g K (e) of Eq. ( 13~TT) . which we write g(e) in this section. 

By designating the Strutinsky smoothing procedure of Eq. ffl5|) with S (do not confuse 
with the S'-matrix that does not appear in the followings), we write g(e) as S[g](e). We 
do not write Sf^e)] since S is not a function but a functional. Using this S operation, we 
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define the result of the application of the reference density Strutinsky method to g(e) by 



Owing to the linearity of the Strutinsky smoothing procedure, the right-hand side of 
Eq. fl53|) can be rewritten as 



which means that SW^] and S[g] differ only where S[g re f](e) ^ fiwte), i-e., where g r ef(e) 
cannot be approximated very well by a polynomial of order 2p over an interval of a few 7 
width. If one defines g re f by oversmoothing g (e.g., g TC f = S[g] with large 7), the reference 
density method and the original Strutinsky smoothing give very close results (i.e., almost 
SVeffg] = S[g}). If one superimposes the peak at energy zero to this reference density, one 
will obtain SW^] which is almost equal to g Te { near energy zero and is very close to S[g] 
anywhere else. 

It should be mentioned that the new smoothing procedure is generally more time- 
consuming than the original one because the integration of g TC f(e) with respect to the 
single-particle energy necessary to calculate S[g re f](e) cannot be done analytically in gen- 
eral. In practice, we sample g re f(e) at an interval of 0.15hu and use a polynomial inter- 
polation between the sampling points. 

B. Construction of the reference density 

The remaining problem is how to determine the shape of the peak at energy zero, 
which is the only important part of the reference density g re f(e). First we present our best 
method. Second, we discuss shortcomings of some other methods which we have tried. 

The best method is to define the reference level density as 



This is the same as the Strutinsky smoothing without the curvature correction polynomial 
(p = 0) except that 7 is a function of energy e, for which we assume 



S re {\g)(e) =S\g- 0rcf](<O + fltef(e)- 



(53) 



StMif) = %]( e ) - %rcf](e) + 0ref(e), 



(54) 




(55) 




(56) 
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FIG. 17: 7(e) defined by Eq. to be used to construct the reference density. 

with parameters 7 a = 3, 7b = 0.7, 7 C = 2, and e pea k = in units of Too. Eq. (!56|) is shown 
graphically in Fig. [T7] to elucidate the roles of each parameter. 

The quantity 7 a is chosen to be large enough so that it holds S^gw] — g r ei and thus 
<Sref[<?] = S[g] at energies distant from e pea k. The quantity 7b is chosen to be small enough 
so that it holds S^ef] — S[g] and thus S^fM = <7ref at energies near e pea k, but not 
so small as the peak is split into more than two peaks. The quantity 7 C is determined 
empirically to obtain smooth results. The energy e pea k is taken as zero for neutrons, while 
it should be around the Coulomb-barrier-top energy for protons. In this paper, we employ 
the reference density method only to treat the neutron spectrum, since for protons the 
standard Strutinsky method works rather well due to the fact that the peak energy is 
considerably larger than the proton Fermi energy even near the proton drip line. 

In Fig. [TU smoothed level densities are shown for the neutron spectrum of 166 Er. The 
result with the reference density method (S Ie { [g], long-dash line) has slightly sharper peak 
at around —1 MeV than the result without using it (S[g], solid line). Their difference 
Ag = S rc f[g] — S[g] looks small but turns out to play an important role in improving 
the plateau in Sec. IIII CI The difference multiplied by ten (10 x Ag) is shown with a 
short-dash line. From Eq. (I54p . Ag = g rc { — 5[g re f], where g Te { is shown with a dot line. 
From the shape of Ag one can see that the Strutinsky smoothing smears the peak of g ve { 
at —1 MeV by moving the density near the top to its hillsides around —7 MeV and 4 
MeV and that the reference density method cancels out this movement by using this Ag 
correction term. 

Incidentally, if one makes 7 a function of e, not of e', in Eq. (IB3|) . one finds fake dips in 
both sides of the peak, as well as an enhancement of the peak. (Changing 7 as a function 
of e in the ordinary Strutinsky method also leads to similar fake dips and bumps.) Our 
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FIG. 18: The Kruppa's level density for a neutron in the ground state of 166 Er smoothed in 
various ways. The basis is specified by = 20. Deformation parameters are /?2 = 0.280 

and /?4 = 0.005. In the Strutinsky smoothing, 7 = 1.2hu and p = 3 are used. See text for 
explanations. 

choice does not suffer from this problem. One should also note that, although the total 
number of levels of the reference level density is not exactly equal to that of the original 
discrete spectrum, 

/oo poo 
g Ie f(e)de ^ / g(e)de, (57) 
-00 J —00 

it still holds 

/oo roc 
S iei [g](e)de = g(e)de. (58) 
-00 J— 00 

We have also examined the possibility of least-square fittings of trial functions. For 
deformed nuclei, we could determine clearly the center and the width of the peak. For 
spherical nuclei, however, we could not because the levels are multiple degenerated and 
thus are very sparse. 

An alternative method to determine the reference density g Te f is the semiclassical es- 
timation. We have tried the level density obtained with the OBTF, only to find much 
stronger dependence on 7 than that of the standard Strutinsky method for a test calcula- 
tion without spin-orbit force. The OBTF approximation does not seem to be sufficiently 
precise to calculate the shell correction energy. Indeed, it is known that one should include 



higher order approximation t 



lan the Thomas-Fermi approximation in the semiclassical 



Wigner-Kirkwood expansion 38|, l39|] . For this purpose, however, one must extend it to 
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the case of truncated oscillator basis expansion. In addition, semiclassical approaches are 
more difficult to use for deformed nuclei. They also have some problems near the particle 
threshold (drip lines) 26]. 



C. Improvement of plateau condition 
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FIG. 19: The neutron shell correction energies for 166 Er and 226 Er calculated with the Kruppa 
method (left), and the Kruppa + reference density (right) Strutinsky method as functions of 
the scaled smoothing width parameter 7 p (|48p in unit of huj with different choices for the order 
of the polynomial p = 3 — 6. As for the basis size N™ x = 30 is used. 

We show how the reference density method improves the plateau condition in Figs. [HO 
and [20] where the Kruppa's prescription is used throughout. Figure [19] depicts the depen- 
dence of the shell correction energies on the scaled smoothing width parameter y p f|48l) 
for the stable and very neutron rich nuclei, 166 Er and 226 Er, considered as examples 
in Sec. Ill Ft Comparing the results with the Kruppa + reference density method (right 
panels) to those with the Kruppa method (left panels), the dependence on the smoothing 
width is weakened remarkably, although it is not completely satisfactory in the case of 
226 Er. The dependence on the order p of the smoothing function is also greatly reduced 
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FIG. 20: Comparison of the neutron shell correction energies as functions of the scaled smooth- 
ing parameter j p calculated with the reference density smoothing (solid curves) and with the 
ordinary Strutinsky smoothing (dashed curves) for 48 Ca, 90 Zr, 146 Gd, 208 Pb, 242 Pb, and the 
superheavy nucleus 298 114, all of which are spherical nuclei (/3 = 0). The order of smoothing 
function is p = 3 and N™ x = 30 is used. 



and the difference of the shell correction energies between p = 3 and 6 is typically within 
a hundred keV. Since this is a general tendency, we show only the results with p = 3 in 
the followings. 

Combined with the reference density method, a better stability against the width 
and the order of the smoothing function is obtained: In most cases, the shell correction 
energy has a minimum as a function of the smoothing width parameter in the range, 
7 P = (l — 2)^a;, and around the minimum we often find a plateau-like quite flat landscape. 
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According to Ref. 25J, the shell correction energy at this minimum should be adopted 
even in the case with no pronounced plateau (the local plateau condition). We show 
examples for several test cases from light to heavy spherical nuclei in Fig. [20j where, in 
each panel, the solid curve is the result with the reference density method and the dash 
curve is the one without it. Comparing two curves, one can see that the plateau is always 
improved by the reference density method. The improvements are remarkable especially 
for 48 Ca and 146 Gd. It should be emphasized that these improvements do not change 
when increasing the size of the model space. However, there are some exceptions as is 
shown for 242 Pb in Fig. [20j where no minimum but a inflection point appears. Although 
the dependence is reduced, it is not enough to obtain plateau-like behavior. Note that 
the parameters of the reference density method are fixed to be the same for all nuclei in 
this paper. There is still some room for their further improvement or optimization. 



D. Kruppa-BCS equation 

Whether the pairing correlation is enhanced or not in nuclei near the neutron drip line 
is still an open question. On one hand, high level density near and above the neutron 
threshold is expected to enhance the pairing [l|, • Thicker neutron skin is also likely 
to make the pairing interaction stronger. On the other hand, the radial expansion of the 
single-particle wavefunctions near the Fermi level will weaken the pair-scattering matrix 
elements. There can be a competition between a spatially expanded normal state and 



a com pac t super state 
effect jsj). 
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631 ] (the latter is an manifestation of the pairing anti-halo 



To take into account all of the above effects, one needs to emp 



oy at least mean field 



models in the Hartree-Fock-Bogoliubov (HFB) formalism 4]J, |65[. To mimic them in 
the shell correction method, one has to extend the method in many aspects by, e.g., 
calculating the pairing matrix elements using the wavefunctions, replacing the BCS gap 
equation with the HFB equation, making the radius and diffuseness parameters (R and 
a) not constants but variables to be optimized like deformation parameters (/^ and in 
this paper), etc. Instead of trying to consider everything, we aim at only one thing, i.e., 
the usage of the Kruppa's level density in the BCS calculation. We call this method (or 
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the resulting equation) the Kruppa-BCS method (equation). 

Near the critical point of the transition between the normal and superfluid phases, it is 
necessary to go beyond the mean field treatment by, e.g., the number projection (4] and/or 
its approximate version, the so-called Lipkin-Nogami method or the random phase 

approximation (RPA) method [69]. Generalizations of the Kruppa prescription to such 
treatments are a quite interesting subject. We restrict, however, to the simplest BCS 
treatment in the present work. 

Although some generalizations arc possible, e.g., to the state-dependent pairing inter- 
action 70|], we consider the simplest seniority-type pairing force for the BCS calculation 
in this paper. The following matrix elements are assumed for the pairing interaction Vp a i r , 

(n'\V paiT \ff) = -GS iV 8 jr J\(< ,)/,.(</). (59) 

In the left-hand side, k (for k = represents the label for the time reversal partner of 
the kth eigenstate of the full or the free single-particle Hamiltonian. It holds that e& = e^. 
When k is a label for a free-particle state, e& should be read as e° k . Scatterings from a pair 
of full-Hamiltonian states into a pair of free-Hamiltonian states, and the reverse processes, 
also appear in the Kruppa-BCS equation to be presented later. In the right-hand side, 
G is a constant while / c (e) is a cutoff factor [71], for which we use a different form from 



that of Ref. 



7l| 



1 + erf 



a + a; 



1/2 



1 + erf 



-e + A + A U N 

int 



1/2 



(60) 



where the error function is defined by erf(x) = —= / e~ l dt. We use the cutoff param- 

V 71 " J° 

eters of the pairing model space, A u = Aj = 1.2 hu> and d cut = 0.2 hu. A is the smoothed 
Fermi level defined by Eq. ( IT2"]) . Incidentally, if one uses Abcs (to be defined in Eqs. ( 16"4"|) 
and (165]) ) instead of A in Eq. ( 160]) . one sometimes encounters an instability caused by a 
positive feedback from Abcs, a part of the solution, to the equation through f c (e). 

The energy of the BCS model for a separable interaction like Eq. (IB^]) can be expressed 
as fl, 



E 



BCS 



where the pairing gap A is given by, 

A = ^ 
2 



A 2 

ev 2 (e)g(e)de- — 



f c (e)u(e)v(e)g{e)de. 



(61) 



(62) 
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(The pairing gap for the state % is state-dependent, / c (ej) A, because of the cutoff function.) 
In Eq. (1ST]) and in the followings, the exchange contribution of the pairing interaction to 
the particle-hole channel is neglected. The constraint on the expectation value of the 
number of particles is expressed as, 



oo 

2 



N = / v\e)g{e)de. (63) 

J — oo 

In the above formulation using integrals, one regards the BCS u and v factors as continuous 
functions of the single-particle energy. It should be reminded that these equations are 
not well defined due to the divergence of the level density if the continuum states are 
included. One has to replace the level density g(e) by that of Kruppa in Sec. Ill Dt Thus, 
we naturally define the Kruppa-BCS model as a model obtained by replacing the ordinary 
level density to the Kruppa's one f l3T|) . g(e) =>- g K (e), in Eqs. (ISTf to (163]) . 

In fact the use of the Kruppa level density makes the gap equation ( 162]) convergent 
without the energy cutoff function f c (e). This is because the integral diverges logarith- 
mically as e — > oo if the level density is constant, while it is shown in Sec. Ill El that 
the Kruppa level density g K (e) oc e -1 / 2 , see Eq. (145]) . However, the convergence is slow 
and it is dangerous to rely on it, so that we use the cutoff function ( 160]) in the following 
calculations. 

Considering that only values at discrete points (e = and e°) contribute in the Kruppa 
prescription, the equation to be satisfied by the minimum-energy state is just the standard 

n 

gap equation |4| and the constraint on the number, but with the additional (negative) 
contributions from the free spectra: 



i M 

N = - V 
2& 



fM) 2 /c(c? 



0\2 



ei - Abcs) 2 + /c(e,) 2 A 2 - A BC s) 2 + /c(e°) 2 A 2 

£j — Abcs e° — A B cs 



(64) 
(65) 



(d - Abcs) 2 + /c(e,) 2 A 2 y/@ - A BC s) 2 + U4) 2 ^ 2 . 

As in the case of the usual BCS equation, the pairing gap and the chemical potential 
(A, Abcs) are determined by these two coupled equations for given force strength G. Note 
that M is two times the number of the pairs of time-reversal states, namely the degeneracy 
is explicitly counted in the level density. One should assume €\ = e 2 < e 3 = e 4 < • • ■ < 
€m-i = &m and e? = e° — e 3 = e 4 — ' ' ' — e °M-i = € °m t° understand the reason of 
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appearing a factor | in many of the equations in this paper, e.g., in the right hand side 
of Eqs. fl64|) and fl65|) . In the case of odd particle number the blocking BCS calculation 
should be done |4|; i.e., the single-particle level occupied by the last odd particle, for 
example i = N, should be eliminated from the pairing model space, and the resultant 
BCS equation with number N — 1 is the same as in the case of even particle number. 

It is known that the BCS equation does not necessarily have finite pairing gap solutions 
if €n < €n + i |4(; namely the system is not in the superfluid phase but in the normal phase. 
The critical force strength G cri t (A = if G < G crit ) is given by 



mm 



G C rit e N<^'<< E N+i 



i f ( fM) 2 



2 £ VI*- -VI k°-A' 



(66) 



where the minimum value of the right hand side is searched with respect to A' (in the 
case of odd particle number, it is always = ejv+i, and for the blocked level i = N 
the minimum value in Eq. (166]) should be searched for e^v-i < A' < ejv+i)- Although the 
Fermi energy A in the normal phase is arbitrary within e^v < A < ejv+i, it is desirable 
to define the Fermi energy uniquely in the later discussion (see Sec JIII H| . Therefore, we 
define Abcs when A = as the A' that gives G CT i t in Eq.(l66lh 

In selfconsistent methods, one only needs to deal with particle-bound nuclei with neg- 
ative Fermi energies. In shell correction approaches, however, one needs some reasonable 
solution for positive energy Fermi levels. This is because negative Fermi levels of the mi- 
croscopic part do not always mean positive separation energies calculated from the total 
energies of the shell correction method, see Sec. IIII HI 

One must be careful in applying the Kruppa-BCS method to the particle-unbound 
cases. More precisely, if the Fermi energy is higher than the lowest energy of the free 
spectra {e°; i = 1, M}. For example, if < < e^+i the right hand side of Eq. ( ]6"S|) 
has no minimum because of the negative contribution of the free spectra, and G cr it cannot 
be defined. As an another peculiar feature of the Kruppa-BCS equation, the solution is 
not always unique for a positive Fermi level. This nonuniqueness is easy to explain for 
the normal states (A = 0). There are more than one ways to fill the spectrum as normal 
states, i.e., to choose the Fermi level A such that ej < A < e^+i, e° < A < e° +1 , and 
i — j = N. For example, in a case e° < < eg < cat +6 , both A = cat+4 and A = tN+6 
have correct number of particles. One has to pay attention to choose the physically most 
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reasonable solution; e.g., the one which gives more continuous total energy with respect 
to the change of deformation parameters. 



E. Extension of the Kruppa's prescription to other observables 

Subtraction of the free contributions in the Kruppa level density in Eq. fl3T|) reminds 
us of the counter term in the renormalization procedure; both contributions diverge but 
the difference remains finite being independent of the cutoff. Therefore, it may be natural 
to extend this idea to other observables: 

(O) =► (0) K = (O) - (O) , (67) 

where the first term is the expectation value with respect to the wave function calculated 
by the diagonalization of the Woods-Saxon potential and the second term is that of the free 
Hamiltonian (or the repulsive Coulomb Hamiltonian for protons). We consider only one- 
body observables in the mean field approximation for the many-body wave function. In 
the simple independent particle approximation, e.g., the Hartree-Fock theory, the second 
(free) term does not contribute as long as the Fermi energy is below the particle threshold. 
If the residual interaction is included, however, the occupation probabilities of unoccupied 
states become non-zero, and then the free-spectrum terms do contribute, which is exactly 
the situation in the case of the BCS theory for the pairing correlation. 

The Kruppa-BCS gap and the number equations, Eqs. (|64j) and (|65|) . can be regarded 
as examples of the above extended procedure ( 1671) because they are derived from 

A = G(P f ) and N = (N), (68) 

where P^ is the pair transfer operator, whose matrix elements are (ii'|Pt|0) = 5m'/ c (q), 
and N is the nucleon number operator. 

Note that the simple BCS calculation of observables composed of the spatial coordi- 
nate r diverges as the basis size is increased for nuclei near the particle threshold. This 
is because the continuum states have finite occupation probabilities; the so-called the 



"neutron gas" problem 4l|, |65| . Therefore, it is impossible to obtain a reliable estimate 



for, e.g., the root mean square radii 72] or the quadrupole moments. It can be shown, 



however, that with the prescription fl67j) such observables also converge as iV™^* — > oo, by 
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employing the oscillator-basis Thomas-Fermi approximation in Sec. Ill El as in the same 
way as for the level density Thus, the Kruppa method relieves the conventional BCS 
method from the failure of the neutron gas problem in nuclei near the drip line. The 
results will be reported elsewhere [3]. 



F. Determination of the strength of the pairing interaction 



The strength of the pairing interaction G is often determined so as to reproduce the 
empirical smooth trend of the pairing gap in the continuous spectrum approximation, 
in which the smoothed level density is used in the BCS calculation js), [lo|, 74-76]. This 
method is almost indispensable in order to treat, say, all the nuclei in the nuclear chart 
on a single footing. A consistent usage of the Kruppa's level density also applies to this 
procedure. 



However, in most of the existing shell correction calculations, e.g., Ref. la, |77J], this 
procedure is not followed rigorously: The so-called uniform level density approximation j^, 
lOf ] is additionally employed, i.e., the energy dependence of the level density is neglected 
and it is replaced by a constant value at the Fermi energy, (/(A). 

As is discussed in Sec. HID} the Kruppa's level density has a peak near the particle 
threshold. We have found that this peak strongly affects the pairing correlation in nuclei 
near the drip line. Therefore, the usual method to approximate the level density as a 
constant, g(X), over the entire energy interval where the pairing is active is inadequate. 
Instead, the energy dependence of the level density should be evaluated exactly. We solve 
the following continuous version of the gap equation and the constraint on the number, 
2 1 r°° f c (e) 2 



G 
N 



1 r - x 
2 



>-A B cs) 2 + /c(e) 2 A 2 
e — Abcs 



g(e)de. 



(69) 



(70) 



'(e - Abcs) 2 + / c (e) 2 A 2 . 
Substituting A with a value from some empirical formula for the pairing gap, one can 
determine the Fermi level Abcs from Eq. (ITDl) and then the force strength G from Eq. (1691) , 
and obtain the smoothed BCS energy, 



E 



BCS 



A 



BCS 



e-A B cs) 2 + /c(e) 2 A2j 



eg{e)de 



A 2 
G ' 



(71) 
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This completes the formula for the total energy in Eq. (|T4|) . The force strength G de- 
termined by Eq. f )69|) is used in the Kruppa-BCS (or the usual BCS) method in the 
previous section. Needless to say, the level density should be replaced, g(e) =>- g K (e), in 
Eqs. ( 1691) to ( I7TT) for the Kruppa-BCS calculation. Although various kinds of input A can 
be presumed 77| . We use a standard choice A = 13/yA MeV in this paper. 

In the uniform level density approximation, the energy integral in Eqs. ( 1691) to ( ITTj) 
can be performed analytically (4] with a sharp cutoff of a pairing model space, A — A] < 
e < A + A u , (d cut — > in Eq. (!60|) ). and the pairing strength can be calculated by 

A, 




(72) 



1 ~ fA x A n 
« - A) log — — 
2 yv ; I A 2 

Moreover, the smooth pairing energy is replaced by the corresponding approximate ex- 
pression; 

1 r 



BCS 



E. 



s.p. 



~A 2 ^(A) 




(73) 



1 7 



4 



A^(A). 



In Ref. 



the results of the continuous BCS equation, Eqs. fl69|) to f |7"Tj) . and of its 
uniform level density approximation, Eqs. ( I72|) and (!73|) . were compared. The difference 
in the smooth pairing energy, -Ebcs — -^s.p. ; was found to be smaller than a few hundred 
keV in most cases. In our calculations, the difference is even smaller, less than a hundred 
keV both for neutrons and protons. 

According to the sharp cutoff in the uniform level density approximation, the number 
of levels included in the conventional BCS calculation (correspon ding to Eqs. (|62|) and 
f )63|) ) is often restricted from i = N\ up to N u in the following way 

n ^r + ~r^^^ m f N -g^Ax+l tor N >g{X)A h 
N u = N + g{\)A u , Ni = < 

[l for N < g(X)Ai. 

Namely, the actual cutoff is often done not for the single-particle energy but for the 
number of levels. In contrast, in the Kruppa-BCS method, the cutoff must be done in 
terms of energy. We use the same cutoff parameters as used in the smooth energy cutoff 
factor fl60|) (A[ = A u = 1.2 fiuj) when we test the sharp cutoffs in Sec. IIII Gl 



(74) 
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G. Results of the Kruppa-BCS calculations 



In this subsection we compare the results of several variants of the BCS calculations. 
We make a combined use of four kinds of classifications to specify those variants. The first 
kind of classification is the Kruppa-BCS method or the ordinary BCS method. The sec- 
ond one is the smooth-energy cutoff or the level-number cutoff. Because the level-number 
cutoff is not applicable to the Kruppa-BCS method, there are three possible combinations 
of the first and the second classifications, which we call "Kruppa-BCS" (with energy cut- 
off), "BCS with energy cutoff", and "BCS with number cutoff". The third classification 
is whether one determines the strength G by the continuous gap equation (|69l) or by its 
uniform level density approximation ( 172|) . We call the former "continuous G" and the 
latter "uniform G" . The fourth is whether one uses the smoothed Kruppa level density 
g K {e) or the usual one g(e) in determining G. There are four possible combinations of the 
third and the fourth classifications, which we call "continuous G with g K " , "uniform G 
with g" , etc. In total there are twelve possible variants, among which we choose the most 
reasonable three to show in Fig. [21] and three unconventional variants to show in Fig. [23] 

In Fig. ED we show the calculated strength G and pairing gap A for neutrons in the 
Z = 68 (Er) isotope chain covering a few more numbers beyond the proton and neutron 
drip lines. The deformation parameters 02 an d 04 are determined to minimize the total 
energy for each nucleus. The basis size is changed as A^™^ x =12, 20 and 30. The 
figure includes the results of three variants of the BCS calculation, A: Kruppa-BCS + 
continuous G with g K , B: Kruppa-BCS + uniform G with g K , and C: BCS with number 
cutoff + uniform G with g. The choices A and B are based on the Kruppa's prescription 
and are new variants we introduce in this paper. The choice C is the conventional one 



employed in, e.g., Refs. [15l. \77\. 



The results of the Kruppa-BCS method (A and B in Fig. [21] (a) to (f)) are, first of all, 
stable against the increase of N™ x as the shell correction energy is. In this way, the 
microscopic quantity -E pa ir, as well as E^, can be calculated to any desired accuracy by 
increasing the basis size. The whole procedure is consistent and unambiguous, which is 
the first and most important purpose of the present paper. Second, in Fig. [21] (a) to (c), 
one sees that the continuous strength G (A) is systematically larger than the uniform 
strength G (B) on the neutron-rich side. This difference can be traced back to the behavior 
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FIG. 21: Neutron's pairing force strength G n (left) and the pairing gap A n (right) as functions 
of the neutron number N for Z = 68 (Er) isotope chain. The top, middle, bottom panels are 
for N™ x = 12, 20, 30, respectively. The results with three different variants of BCS calculations 
are included; A: Kruppa-BCS + continuous G with g K , B: Kruppa-BCS + uniform G with g K , 
and C: BCS with number cutoff + uniform G with g (see text for a detailed explanation). The 
solid curves, G = 20/ A [1/MeV] and A n = 13/-/A [MeV], are also included. 

of the Kruppa's level density in Fig. El enlargement of which in the energy range of the 
most influential part to BCS calculations is shown in Fig. [221 The Kruppa level density 
decreases as the single-particle energy exceeds the threshold, which leads to the increase 
of the pairing strength compared to the case of uniform level density for nuclei near the 
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FIG. 22: Enlargement of the smoothed level densities with and without Kruppa prescription 
shown in Fig. [3] in the energy range of the most influential part for BCS calculations. The 
Kruppa's densities calculated with N™ x = 20 and N^ x = 30 are almost indistinguishable. 



drip line (note that the strength is inversely proportional to the level density at the Fermi 
level). This means that the uniform level density approximation is inappropriate when 
the peak at energy zero is close to the Fermi level. Third, in Fig. [21] (d) to (f), one sees 
that the pairing gap of the calculation A is closer to the input value A = 13/y/A MeV 
than that of B, which means that the continuous G choice is preferable to the uniform 
G choice. We propose the method A, i.e., the Kruppa-BCS method with the strength G 
calculated by the continuous gap equation with the Kruppa's level density, as the best 
method for reliable calculations of, e.g., nuclear masses. 

Let us examine also the results of the conventional BCS calculation. When the basis 
is as small as N™ x — 12, the pairing gap of the conventional BCS calculation (C in 
Fig. [21] (d)) agrees very well with that of the Kruppa-BCS with continuous G (A). This 
agreement is, again, accidental since increasing the basis size changes the results consid- 
erably (compare the calculations C in Fig. [21] (d) to (f)). Even when the basis size is 
as large as N™ x — 30, the pairing gap of the conventional treatment (C in Fig. [21] (f)) 



does not look totally wrong 



16], e.g., being closest to the input values for iV < 132. 



However, the force strength G in this case (C in Fig. [21] (c)) behaves rather unnaturally. 
It is quite different from those of A and B in Fig. [21] (c) as well as from frequently used 
simple expressions like G = 20 /A [1/MeV]. It is large in the stable region but decreases 
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FIG. 23: Same as in the right panels of Fig. [21] but with non-conventional choices of the BCS 
calculation; D: BCS with energy cutoff + uniform G with g, E: BCS with number cutoff + 
uniform G with g K , and F: BCS with energy cutoff + uniform G with g K (see text for a detailed 
explanation) . 



dramatically toward the neutron drip line. This peculiar behavior is caused by a spurious 
effect caused by including more and more continuum states coming into the pairing model 
space when increasing the basis size in the conventional treatment with g (not with g K ), 
as is clearly seen in Fig. 
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This strong reduction in G combined with the level number cutoff is helpful to prevent 
the pairing gaps from becoming extremely large as is often the case for the other unrea- 
sonable variants (see Fig. 123]) . However, this result is also unphysical as it is inspected 
from the results of non-conventional calculations shown in Fig. [23] Here we show three 
non-conventional choices of the BCS method, D: BCS with energy cutoff + uniform G 
with g, E: BCS with number cutoff + uniform G with g K , and F: BCS with energy cutoff 
+ uniform G with g K . The cases E and F are included to see what happens if reasonable 
values of the strength G (calculated by g K instead of g) are used. For sufficiently large 
basis sizes, the resultant pairing gaps are too large, if the energy cutoff (D in Fig. [23] (c)) 
is used in place of the level number cutoff (C in Fig. [21] (f)). Other non-conventional 
choices (E and F in Fig. 1231) uses more reasonable pairing force strength G (the same as B 
in the left panels in Fig. |2~T|) . and give reasonable values of A in stable nuclei [N < 120), 
but the pairing gaps diverges when approaching to the drip line (see Fig. [23] (c)). Only 
with relatively small basis sizes such as iV™^ x = 12, reasonable pairing gaps are obtained; 
actually all the six variants, A to F in Figs. [21] and [23] gives almost the same results for 
stable nuclei with iV ™ ax = 12. 

A lesson of these test calculations is that the conventional BCS calculation is very 
dangerous if one uses the continuum states obtained by the diagonalization with a large 
basis. The subtraction procedure of the free contributions, i.e., the Kruppa prescription, 
is indispensable to treat the pairing correlation in nuclei far from the stability. 



H. Readjustment of the potential depth for the Fermi level consistency 

In Strutinsky calculations, the neutron and proton Fermi levels of the single-particle 
potentials are not equal to the derivatives of the total energy, dE/dN or dE/dZ, in 
general unlike in mean field models. As it is explained in Sec JIIBl the total energy is 
divided into the macroscopic and microscopic parts and they are calculated separately: 
The Fermi energies corresponding to them are generally different. Since most of nuclei 
are in the superfluid phase and the microscopic energy is calculated by the BCS method 
(see Eq. ([14]) ). we define the Fermi energies, 

_ dE mac (N,Z) _ <9ff mac (iV,Z) 

K ~ M , A p - — , (75) 
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for the macroscopic part, and for the microscopic part 



\ mic \ 

A„ = A 



(n) 
BCS 



\ mic \ 

A„ = A 



(p) _dE$ s (Z) 

BCS — 



(76) 



ON ' p ~ ^ b dZ ' 
The "total" Fermi energy is directly related to the two-particle separation energy, — 5 < 2 n /2 
or — 5 < 2p/2, with the definition of the separation energies, 



(77) 



S 2n (N, Z) = £(iV - 2, Z) - E(N, Z), 
S 2p (N, Z) ee £(iV, Z - 2) - £(AT Z). 

Although the physically meaningful quantity is only the total one, it is desirable that 
all the macroscopic, microscopic, and total Fermi energies coincide with each other. As 
it is shown in the following, however, most of the existing Woods-Saxon parameter sets 
lead to A™ 1C > and A™ 1C > 0, i.e., particle- unbound, at the drip lines, which is not only 
inconsistent conceptually but also can be problematic for the Kruppa-BCS calculation as 
is mentioned in Sec. IIIIDI Therefore, we consider how to avoid this problem. 
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FIG. 24: Two-nucleon drip lines calculated with the universal parameter set. Solid lines are 
two-nucleon drip lines, which pass between adjacent bound and unbound even-even nuclei. The 
dash lines are the boundary between nuclei having the positive and negative microscopic Fermi 
levels. The dot lines are the drip lines of the macroscopic part of the model. N™ x = 20 is used. 
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As a test calculation of the improved microscopic-macroscopic method developed in the 
present work, we have done global mass calculations for even-even nuclei with 8 < iV < 
184 and 8 < Z < 126. The basis size specified by N™ x — 20, the Strutinsky smoothing 
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FIG. 25: Same as in Fig.[23]but for other parameter sets. Only the neutron drip lines are shown. 
The ordinate values are shifted upward by 30n where n = for the parameter set Wahlborn, 
n = 1 for Rost, n = 2 for Chepurnov, n = 3 for Wyss-1, and n = 4 for Wyss-2. 



parameters 7 = 1.2 hu and p — 3, and the smoothed pairing gap A = 13/ v A MeV are 
used for these mass calculations. In Fig. [2U we show the two-neutron drip line (& r = 0) 
and the two-proton drip line (^p = 0) calculated with the universal parameter set |43l. 144 1 
for the Woods-Saxon potential. The drip lines tend to fluctuate outside the drip lines of 
the macroscopic part of the model, defined by equations A™ ac = or A™ ac = 0, due to the 
shell effect. On the other hand, the line in which neutron's Fermi level is zero (A™ 10 = 0) 
is located by 7 (12) neutrons inside the line A™ ac = at Z=A0 (80): The microscopic 
Fermi energy A™ 10 is positive and non-negligible at the neutron drip line. 

Neutron drip lines for five other potentials are shown in Fig. [25j For potentials of 
Wahlborn [79] and Rost [80[, the displacement of the line S2 n = from the line A™ 1C = 



is as large as AiV=6 to 9 (13) at Z=A0 (80). For potentials Wyss-1 |4(j and Wyss-2 [81] 
AiV=2 (9) at Z=40 (80). (The values of the parameters for the potential Wyss-2 can be 



found in Table I of Ref. 



821 ] . The numberings for Wyss's two potentials are tentative. 
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The smallest displacement is obtained for Chepurnov's potential [83| for which AiV =2 (5) 
at Z=A0 (80). These displacements clearly show that the linear dependence of the depth 
of the potential on iV — Z as in Eq. (j3J) is oversimplified for nuclei far from stability. 

As for protons, on the other hand, the situation is much better; three lines almost 
coincide in Fig. [21] as well as for the other potentials (not shown), which will be mainly 
because the proton drip line is closer to experimentally known nuclei than the neutron 
drip line is. However, again, the microscopic Fermi energy A™ 1C > at some places on the 
proton drip line. 

The reliability of shell correction energies for nuclei in the area between lines A™ 1C = 
and = is not very high because E g _ p _ is affected by basis-dependent discretized con- 
tinuum levels directly (not via smoothing). Thus, it is preferable to modify the potential 
parameters in such a way that the Fermi levels are consistent with the liquid-drop part of 
the model. Among the parameters of the central potential in Eq. ([2]), the radius -Roce and 
the surface diffuseness acE are related directly to other observables than energy. Hence 
we choose the depth, 



N- Z 

1 ± K C E 



A 



(78 



Vdepth — —Voce 
(see Eq. P|) to modify. 

Our procedure is as follows. We consider only spherical shape (we use the same depth 
for deformed shapes), and we neglect the spin-orbit potential. The single-particle Hamil- 
tonian is then given by Eq. ( 1331) . The local number density of neutrons or protons at 
the energy e is represented by p TF (r, e) of Eq. (1351) in the Thomas- Fermi approximation. 
Given a set of neutron and proton numbers (N, Z), we substitute e with A™ ac in Eq. (J75J) 
and readjust Vd ep th so as to fulfill 

N = 4:71 p T Jr, \™ c )r 2 dr. (79) 
Jo 

For protons, N and AJf ac are replaced with Z and A™ ac , respectively, and the integral is 
only inside the Coulomb barrier for A™ ac > 0. It should be noted that, in our method, the 
depths of the central potentials are determined by the other parameters than Voce and 
^ce- We use the original parameter value of Voce (multiplied with Aso) only to determine 
the depths of the spin-orbit potentials. We do not use kqe anywhere. 

If the Fermi level of the liquid-drop model is close to zero, the readjusted depth of the 
potential becomes significantly shallower than the smooth continuation from the results 
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FIG. 26: The depths defined in Eq. (|78p of the central potentials of the universal parameter set 
(solid and dash lines) and those readjusted in the Thomas-Fermi approximation (filled markers) 
for Er isotopes. For nuclei near or beyond the particle thresholds (or the Coulomb barrier top 
for protons) extrapolated values are plotted (empty markers). N™ x = 30 is used. 

for more bound nuclei. It is due to the tail of the Woods-Saxon potential. This deviation 
from the smooth trend does not seem to be physically meaningful. Thus we switch to an 
extrapolation of the smooth trend if the Fermi level of the liquid-drop model is higher 
than some predefined energy: We take this energy as —2 MeV in this paper. For neutrons 
(protons), we use a polynomial of second degree in N (Z) determined by three heaviest 
even-iV isotopes (even-Z isotones) not matching the above condition. 

This extrapolation is also indispensable to determine the potential depth for nuclei 
which are outside the drip lines of the liquid-drop model. One has to calculate such nuclei 
because they may be bound since the shell effect can shift the drip lines. It is also because 
nuclei just beyond the drip lines are necessary to determine the drip lines themselves. 

The resulting readjusted potential depths are shown in Fig. [26] for Er isotopes. The 
original parameter set is the universal one. While Vd ep th is readjusted, the other pa- 
rameters are kept unchanged. For both proton and neutron, the changes due to the 
readjustment are almost vanishing for stable nuclei at iV ~ 100: This is totally non- 
trivial because the original parameter is determined by completely different requirements. 
In the neutron (proton) drip line at N = 156 (N = 76), the readjustment of the neutron 



57 



5 


> 5 

a; -j 

f -10 

■ r— 1 

g -15 
-20 
-25 

60 80 100 120 140 160 

AT 

FIG. 27: Proton and neutron Fermi levels of Er isotopes calculated using the universal parameter 
set with and without the readjustment of the depth of the central potentials (four kinds of lines). 
Two nucleon separation energies divided by two with the sign inverted are also shown (four kinds 
of markers). N™ x = 30 is used. 

(proton) potential has non- negligible size AVdepth = —2.3 (—1.3) MeV. This means that 
the original parameter set is quite reasonable near the /3-st ability line but not sufficiently 
accurate to be applicable to the driplines. 

Our readjusted depths look very smooth functions of N and Z, which seem to be fitted 
nicely with simple functions having only a few parameters. Such a fitting ought to be 
done when we will publish an optimized parameter set in future. Indeed, Nazarewicz et 
al. introduced an extra (N, Z) dependence to potential parameters for the same purpose 
[lij. At present we make a direct use of the depths determined in the Thomas- Fermi 
approximation. 

In Fig. [27] we show the various kinds of Fermi levels for Er isotopes. One can see near 
drip lines that Fermi levels calculated with the readjusted potential depths are in good 
agreement with separation energies. This fact, as well as the fact that the original and 
readjusted depths are almost equal for stable nuclei, confirm the soundness of our pre- 
scription, although it may change slightly the single-particle spectrum from the optimized 
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FIG. 28: Same as in Fig. [23] but with readjusted depths for the central potentials. 
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FIG. 29: Same as in Fig. [25] but with readjusted depths for the central potentials. 



one of the original potentials. 

The drip lines with the readjusted potential depths are shown in Fig. [28] for the uni- 
versal parameter set and in Fig. [29] for the other parameter sets. One can see that the 
three lines are now quite close to each other. With this readjustment method, one can 
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treat more reliably nuclei near the driplines. Furthermore, it is worth stressing that one 
can control the drip lines by changing the macroscopic part while keeping automatically 
the consistency with the microscopic part. 

Incidentally, the relation between macroscopic part and the microscopic part has been 



payed attention to by Myers 42] already in 1970. He determined the droplet model pa- 
rameters in terms of a Thomas- Fermi statistical model with a phenomenological velocity- 
dependent force applied to infinite and semiinfinite nuclear matter. However, this ap- 
proach seems rather distant from what is proposed in the present paper. 

IV. CONCLUSION 

We have examined the Fourier transform of the smoothing function of Strutinsky to 
find that it is nothing but a low-pass filter, which passes only short-time components. The 
polynomial part of the filter is simply a truncated Taylor expansion of e^'/ 2 ) 2 where k is 
the time multiplied by the smoothing width 7 and divided by h. It may be redefined as a 
polynomial to minimize the distortion of the filter near k = 0, which seems simpler than 
the original definition as a curvature correction. We have also derived a relation between 
7 and the order of the polynomial part p that changing 7 proportionally to yjp leaves the 
results of smoothing almost unchanged. This picture of the Strutinsky smoothing as a 
low-pass filter is general and will be useful for investigating the other smoothing functions, 



e.g., those of Ref. [58], than the standard one considered in the present work. 

From this point of view, we have given a negative perspective to the problem of the 
plateau for the Woods-Saxon spectrum, a problem concerning the shell correction energy 
in the microscopic-macroscopic method. It has been known that the shell correction 
energy for the Nilsson spectrum behaves like a long flat plateau as a function of 7 while 
that for the Woods-Saxon spectrum does not show such a magnificent plateau in general. 
We have noticed that the Fourier transform of the Nilsson spectrum has an interval of 
time components where the amplitude is almost vanishing, while that of the Woods-Saxon 
spectrum does not have such an interval. A plateau appears when the cutoff of the filter 
is in such an interval. 

Instead, we have proposed a new method to weaken the dependence on the smoothing 
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width 7. We call it the reference density method, in which the smoothing is applied only 
to the deviation from a reference level density, which was prepared in such a way that the 
peak around energy zero of the Kruppa's level density, which is another principal subject 
of this paper, is not washed away. We have demonstrated that the method works well in 
the desired direction. 

To apply the Woods-Saxon-Strutinsky method (the microscopic-macroscopic method 
with finite-depth potentials) to nuclei near the nucleon drip lines, it seems necessary 
to employ the Kruppa's prescription for positive energy levels, in which the spectrum 
is defined as the Woods-Saxon spectrum subtracted by the free-nucleon spectrum, both 
of which are obtained through diagonalizations in the same oscillator basis. We have 
discussed the ground for this prescription as well as its relation to the continuum level 
density. 

We have also proposed the oscillator-basis Thomas-Fermi approximation, with which 
one can describe spectra obtained from diagonalization in truncated oscillator bases. 
We have demonstrated that this approximation can reproduce average behaviors of the 
Woods-Saxon, free, and Kruppa's level densities. One can also use this approximation 
to show analytically the convergence of the results with Kruppa's prescription versus the 
size of the oscillator basis. 

We have also introduced the Kruppa-BCS method, in which we modified the BCS 
equation for the pairing correlation so that it can be applied to the Kruppa's level den- 
sity by taking into account negative contributions from the free-nucleon spectrum. The 
Kruppa-BCS method is applicable to any cases, while the ordinary BCS method becomes 
very faulty especially when the diagonalization basis is not small and the nucleus is very 
neutron-rich. 

We have also studied how to determine the strength of the pairing interaction to be used 
in the Kruppa-BCS method. An important conclusion is that, in adjusting the strength 
to reproduce the empirical smooth trend of the pairing gap with the smoothed Kruppa's 
level density in the gap equation, one should carry out the energy integral without using 
the uniform level density approximation. 

The inconsistency between the macroscopic and the microscopic parts (i.e., the liquid 
drop model and the single-particle potential) is another problem to the application of 
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the method to nuclei near the drip lines. Calculating masses in the whole nuclear chart 
with several parameter sets for the Woods-Saxon potential using the methods developed 
in this paper, we have found that the neutron drip line of the microscopic part is located 
typically more than ten neutrons inside the dripline of the total energy. We have proposed 
a method to readjust the depths of the central potentials to achieve the consistency within 
the Thomas-Fermi approximation. Although the method contains two simplifications, 
assuming the spherical symmetry and neglecting the spin-orbit potential, the method has 
worked very well to shift the microscopic dripline close to the total dripline. 
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We are going to apply the methods presented in this paper to extend our studies 
on the origin of the prolate-dominance of the atomic nucleus from the Nilsson potential 
to the Woods-Saxon potential. 
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